BEI and equity paper analyses | FINAL | Hurdle models

# Clean up dataframe to keep only meaningful vars and drop NA, ready for LM
tidy_df <- function(src_df, ct, f, drop_islands=TRUE) {
  clean_bei_df <- src_df %>%
    select(CT_UID, all_of(all.vars(f))) %>%
    drop_na() %>%
    units::drop_units()
  
  # Join geom from CT dataframe
  df <- ct %>%
    transmute(CT_UID = GeoUID) %>%
    inner_join(clean_bei_df, by="CT_UID")
  
  # Build neighborhood (with islands) and extract island IDs
  ct_nb <- poly2nb(df)
  
  if (drop_islands) {
    islands <- lapply(ct_nb, min) %>% lapply(function(e) e == 0) %>% unlist %>% which()

    # Cleanup subset DF to drop islands and NA
    df <- df %>%
      filter(!row_number() %in% islands) %>%
      mutate(ct_no = row_number()) # Add CT number, matching row number
  
    # Recompute neighborhood
    ct_nb <- poly2nb(df)
    
    # Clean up final dataframe (no island, no NA's, pure dataframe)
    clean_df <- df %>%
      as.data.frame() %>%
      select(CT_UID, ct_no, all_of(all.vars(f)))
  }

  nbw = nb2listw(ct_nb, zero.policy = !drop_islands)
  nbmatx = nb2mat(ct_nb, style = "B", zero.policy = !drop_islands)
  
  list(df=clean_df, nbmatx=nbmatx, nbw=nbw, island_dropped=drop_islands)
}

This document builds on BEI_equity_paper.Rmd doc, which started as an analysis for the equity paper and moved to an exploration of the best modeling framework (simple OLS, spatial auto-regressive linear regression, spatial random effect model and INLA) to assess the relation between interventions and equity metrics. The conclusions of this exploration led to the following decisions:

1. go with quintiles instead of absolute numbers (easier to interpret)
2. use [zero inflated] Poisson distribution where fit (ie Bike lane ratio)
3. use spaMM package to model association; INLA could be the focus of another methodological paper which would compare both approaches
4. use non-buffered dependent variable, as spatial dependency is already taken care of through the spatial random effect

This document contains the analyses to be included in the BEI and equity paper.

1 Proposal

1.1 Original grant proposal for Census tract level analyses

1.1.1 Objective 1: Socio-spatial inequalities in urban interventions:

• To test if urban interventions tend to be located in low SES neighborhoods (H1a); we will run the following model: \[Log(𝔼(UI_{𝑖} \mid X_{𝑖})) = \alpha + \beta_{1} * SES_{𝑖}\] a Poisson regression model with \(UI\) = the number of Urban Interventions between 2016 and 2021, \(SES\) = the socio-economic status of the census track at 2016.
• To test the reduction in socio-economic inequalities in urban conditions (e.g. length and type of bicycle lanes, NDVI) between 2006 and 2020 (H1b), we will run the following model: \[𝔼(UC_{ij} \mid X_{ij}) = \beta_{0j} + \beta_{1j} ∗ SES + \beta_{2j} ∗ Time + \beta_{3j} ∗ SES ∗ Time + \epsilon_{ij}\] a multilevel linear model with \(UC\) = the Urban condition, \(SES\) = the socio-economic status of the census track.

1.1.2 Objective 2: Examining causal pathways and directionality of intervention/gentrification relation:

• To test if urban interventions occur before gentrification or if gentrification occurs before the implementation of new urban interventions, we will run two models : (1) \[Logit(𝔼(G_{i} \mid X_{i})) = \alpha + \beta_{1} ∗ UI_{i} + + \beta_{2} ∗ PG_{i} + \beta ∗ X_{control_{i}}\] a logistic regression model among low SES census tracts (2006) with \(G\) = Gentrification between 2006 and 2016, \(UI\) = the number of Urban Interventions between 2001 and 2006, and \(PG\) = SES status 2006; and (2) \[Logit(𝔼(UI_{i} \mid X_{i})) = \alpha + \beta_{1} ∗ G_{i} + \beta_{2} ∗ UC_{i} + \beta ∗ X_{control_{i}}\] a Poisson regression model, with \(UI\) = the number of Urban Interventions between 2016 and 2023, \(G\) = Gentrification between 2006 and 2016, and \(UC\) = Urban Conditions in 2016.

1.2 Paper objective

We focus on obj #1: are urban interventions tend to be located in low SES neighborhoods. Here, the urban interventions considered are bike lanes and canopy/tree coverage and low SES ~ high Pampalon deprivation index. In a second step, we will look at the variations of UI and SES and their association.

Data extraction and pre-analyses for the paper looking at BEI and equity. BEI comprise bike lanes and canopy changes while equity is measured through Pampalon deprivation index. (Partially adapted from original work on BEI bike lanes – see bike_lane_stats.R and ReadMe.md.)

Paper is available here.

General processing steps:

• Get CT 2016 boundaries (from Cancensus API)
• Get BEI interventions (canopy / bike lanes) for the selected years (2011 / 2015 / 2016 / 2019)
• Compute changes between years
  • for the original CT boundaries
  • for buffered CTs (250 / 500 / 750m)

In a second phase, these BEI changes will be linked to Pampalon index for 2016 (and 2011 ?)

1.3 Revision history

UPDATE 2021-12-02 Following discussion with @Yan, add normalized bike line changes:

• by CT/buffer area
• by street length within CT/buffer

UPDATE 2021-12-08 Following discussion with @Yan

• Keep only normalized variation of bike lane length
• Reverse relation: BEI metric = f(Equity metric (Pampalon / gentrification))
• Display association at baseline, then add delta over time
• Canopy : distinguish between high and low canopy

UPDATE 2022-01-14 Following discussion with @Ruben and @Yan

• Test model with UC 2011 when modeling UI
• Models with random slope probably not pertinent (see model 6, section 7.1.6)
• Try to define Year variable as continuous instead of category in LME
• Check if range of wSCOREMAT matches theoretical range of x-axis in graph (-.2 >> .2), might ponder the magnitude of the observed trend
• Test three-way interaction instead of 2 times two-way interactions (interactions might be complicated to interpret…)
• For gentrification and multilevel model, use GEE (logit)

UPDATE 2022-02-10 Following gentrification meeting

• Add correlation matrix of SES variables (Caislin)
• Drop wSCORESOC, which does not capture equity dimension (Meghan) [replaced by visible minority]
• Paper focus should be on gentrification

UPDATE 2022-02-15

• Add spatial autocorrelation test (Moran’s I)
• Add Lagrange Multipliers to test best spatial model
• Add comparison of OLS/SAR/SEM models

UPDATE 2022-06-28

• Add categorical SES variables (instead of ordinal coding) to better explain the Q5 vs other effect (see Behzad and Yan’s comments); this is particularly true for bike lanes changes, as no significant effect found yet graph shows clearly a Q5 vs. other effect
• Build maps of SES (?) and BEI / BEC (might be better in ArcGIS)

1.4 A few references…

• General resources on GLMM
  • GLMM FAQ
• Count and rate models:
  • A tutorial on count regression and zero-altered count models for longitudinal substance use data
  • Regression Analyses of Counts and Rates: Poisson, Overdispersed Poisson, and Negative Binomial Models
  • The Analysis of Count Data: A Gentle Introduction to Poisson Regression and Its Alternatives
• Spatial models:
  • https://rspatial.org/raster/analysis/7-spregression.html
  • https://rpubs.com/quarcs-lab/tutorial-spatial-regression
  • https://www.r-bloggers.com/2019/09/spatial-regression-in-r-part-1-spamm-vs-glmmtmb/
  • https://crd230.github.io/lab8.html
  • http://labs.bio.unc.edu/Buckley/documents/AnselinIntroSpatRegres.pdf
• Interactions:
  • Decomposing, Probing, and Plotting Interactions in R
  • Sophisticated models in emmeans
  • R for Researchers: An Introduction
  • Recommendations for presenting analyses of effect modification and interaction
  • Testing interact terms

2 Built Environment Intervention Extraction

2.1 Bike lane changes

We use data categorized by Philippe Apparicio’s team who manually identified bike lanes for each census year since 1991. For this study, we limit ourselves to 2016, 2011 and 2006 census years.

On top of the original CT boundaries, three levels of buffer have been applied to the CT – 250m, 500m & 750m. Then the same series of processing steps (see above) have been applied to the buffers.

# Bike lanes, from Ph. Apparicio
reseau <- st_read(dsn="data/ReseauCyclableFinal.gdb", layer = "Reseau") # Already in NAD83 / MTM zone 8

## Reading layer `Reseau' from data source 
##   `/Users/benoit/WORKSPACE/gentrification_BEI_equity/data/ReseauCyclableFinal.gdb' 
##   using driver `OpenFileGDB'
## Simple feature collection with 82166 features and 72 fields
## Geometry type: GEOMETRY
## Dimension:     XYZ, XYZM
## Bounding box:  xmin: 266985.5 ymin: 5029251 xmax: 320986.1 ymax: 5062652
## z_range:       zmin: 0 zmax: 43
## m_range:       mmin: 0 mmax: 43
## Projected CRS: NAD83 / MTM zone 8

bike_lane <- reseau %>%
  filter(An2016 == 1 | An2011 == 1) %>%
  select(IdRte, ClsRte, Zone, starts_with("An"), starts_with("Typo_")) %>%
  st_cast("MULTILINESTRING") # Get rid of a few MULTICURVE geometries

# CT boundaries for Montreal
CT16 <- get_census(dataset='CA16', regions=list(CMA='24462'), level='CT', geo_format = "sf") %>%
  filter(Type == "CT") %>%
  mutate(interact_aoi = CD_UID %in% c(2466, 2465, 2458)) %>% # Flag Montréal island, Laval and the South shore (Longueuil, St-Lambert, Brossard)
  st_transform(st_crs(bike_lane))

## Reading geo data from local cache.

CT11 <- get_census(dataset='CA11', regions=list(CMA='24462'), level='CT', geo_format = "sf") %>%
  filter(Type == "CT") %>%
  st_transform(st_crs(bike_lane))

## Reading geo data from local cache.

compute_bikelane_by_area <- function(sf_areas, year_fld, typo_fld) {
  # Compute length of bike lanes within each area.
  # --
  # Parameters:
  #   - sf_areas: sf class object defining the areas of interest, must have a GeoUID field
  #   - year_fld: field name specifying the year of interest, e.g. An2016
  #   - typo_fld: field name specifying the typology for the year of interest, e.g. Typo_2016
  
  year_fld <- enquo(year_fld)
  typo_fld <- enquo(typo_fld)
  
  # Compute intersection of bike lanes with areas
  bk <- bike_lane %>%
    filter(!!year_fld == 1) %>%
    st_intersection(sf_areas) %>%
    mutate(bike_lane_length = st_length(.)) %>%
    as.data.frame() %>%
    group_by(GeoUID, !!typo_fld) %>%
    summarise(bike_lane_length = sum(bike_lane_length)) %>%
    ungroup() %>%
    pivot_wider(names_from = !!typo_fld, names_prefix = "Bike_class", names_sort = TRUE,
                values_from = bike_lane_length, values_fill = units::set_units(0, m)) %>%
    mutate(Bike_lane_total = units::set_units(rowSums(select(., starts_with("Bike_class"))), m))
  
  # Merge back into original sf_areas
  bk <- sf_areas %>%
    left_join(bk)
  
  # Replace NA by 0, which occur in Bike_class length
  bk[is.na(bk)] <- 0
  
  return(bk)
}

compute_streetlength_by_area <- function(sf_areas) {
  # Compute length of streets within each area.
  # --
  # Parameters:
  #   - sf_areas: sf class object defining the areas of interest, must have a GeoUID field

  # Compute intersection of streets with areas
  bk <- reseau  %>%
    st_cast("MULTILINESTRING") %>% # Get rid of a few MULTICURVE geometries
    st_intersection(sf_areas) %>%
    mutate(street_length = st_length(.)) %>%
    as.data.frame() %>%
    group_by(GeoUID) %>%
    summarise(street_length = sum(street_length)) %>%
    ungroup()
  
  # Merge back into original sf_areas
  bk <- sf_areas %>%
    mutate(shape_area_km2 = units::set_units(st_area(.), 'km^2')) %>%
    left_join(bk)
  
  # Replace NA by 0
  bk[is.na(bk)] <- 0
  
  return(bk)
}


# Compute year 2016 and year 2011 bike lanes within 2016 CTs
# NB: contrary to the original work, we keep the same area of reference, i.e. 2016
bike_lane_by_CT16 <- compute_bikelane_by_area(CT16, An2016, Typo_2016)
bike_lane_by_CT11 <- compute_bikelane_by_area(CT16, An2011, Typo_2011)
bike_lane_by_CT06 <- compute_bikelane_by_area(CT16, An2006, Typo_2006)

# Compute total street length within CT/buffer
street_length_by_CT16 <- compute_streetlength_by_area(CT16)

# Reorganize data to have all data in one dataframe
bike_lane_changes <- CT16 %>%
  left_join(select(as.data.frame(street_length_by_CT16), GeoUID, street_length, shape_area_km2), by="GeoUID") %>%
  left_join(select(as.data.frame(bike_lane_by_CT16), GeoUID, starts_with("Bike_")), by="GeoUID") %>%
  left_join(select(as.data.frame(bike_lane_by_CT11), GeoUID, starts_with("Bike_")), by="GeoUID", suffix=c(".2016ct", ".2011ct")) %>%
  left_join(select(as.data.frame(bike_lane_by_CT06), GeoUID, starts_with("Bike_")), by="GeoUID")

# Compute ratio of bike lane vs street length
bike_lane_changes <- bike_lane_changes %>%
  transmute(GeoUID = GeoUID,
            interact_aoi = interact_aoi,
            street_length = street_length,
            street_length.hm = street_length / 100, # in hectometers, can be kept as float (not a dependent variable following a Poisson distribution)
            Bike_lane_total.2016ct = Bike_lane_total.2016ct,
            Bike_lane_total.2011ct = Bike_lane_total.2011ct,
            Bike_lane_total.2006ct = Bike_lane_total,
            Bike_lane_total.hm.2016ct = as.integer(round(Bike_lane_total.2016ct/100)),
            Bike_lane_total.hm.2011ct = as.integer(round(Bike_lane_total.2011ct/100)),
            Bike_lane_total.hm.2006ct = as.integer(round(Bike_lane_total/100)),
            Bike_lane.by.street.2016ct = 100 * Bike_lane_total.2016ct / street_length, # In percents
            Bike_lane.by.street.2011ct = 100 * Bike_lane_total.2011ct / street_length,
            Bike_lane.by.street.2006ct = 100 * Bike_lane_total / street_length)

# Compute change between 2011 and 2016 (only for total bike lane length)
bike_lane_changes <- bike_lane_changes %>%
  mutate(Bike_lane_diff.2011.2016ct = Bike_lane_total.2016ct - Bike_lane_total.2011ct)

# Normalize bike lane change by (i) street length and (ii) area
bike_lane_changes <- bike_lane_changes %>%
  mutate(Bike_lane_diff.by.street.2011.2016ct = 100 * Bike_lane_diff.2011.2016ct / street_length)

# Save results
st_write(bike_lane_changes, dsn = "data/bike_length_changes.gpkg", delete_layer = TRUE)

## Deleting layer `bike_length_changes' using driver `GPKG'
## Writing layer `bike_length_changes' to data source 
##   `data/bike_length_changes.gpkg' using driver `GPKG'
## Writing 970 features with 15 fields and geometry type Multi Polygon.

# Clean up
rm(buf_CT16)
rm(bike_lane_by_CT11, bike_lane_by_CT16, bike_lane_by_CT06)
rm(street_length_by_CT16)

2.1.1 Illustration of bike lane metrics

Check output for one specific dataset (Census tracts 2016, no buffer)

2.1.1.1 Raw length

Bike lane length in 2016 within CTs, measured in meters

ggplot() +
   geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_total.2016ct)), lwd=0) +
  scale_fill_continuous(name = "Total length (m)")+ 
  labs(title = "Length of bike lanes within 2016 CTs")

2.1.1.2 Absolute length change

Absolute bike lane length change between 2011 and 2016, in meters

ggplot() +
   geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_diff.2011.2016ct)), lwd=0) +
  scale_fill_gradient2(name = "Length (m)")+ 
  labs(title = "Changes in bike lane between 2011 and 2016")

2.1.1.3 Normalized length change (by street)

Relative bike lane length change between 2011 and 2016, normalized by street length within CT in 2016, expressed in %

\[Variation = \frac{Bike Lane_{2016}}{Street Length_{2016}} - \frac{Bike Lane_{2011}}{Street Length_{2016}}\]

ggplot() +
   geom_sf(data=filter(bike_lane_changes, interact_aoi), mapping = aes(fill=as.numeric(Bike_lane_diff.by.street.2011.2016ct)), lwd=0) +
  scale_fill_gradient2(name = "Variation (%)")+ 
  labs(title = "Changes in bike lane between 2011 and 2016, normalized by street")

2.1.2 Basic stats on each layer

At the Census Tract level.

# CT level
ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.2011.2016ct)), binwidth = 250) + 
  xlab("Difference of bike lane length between 2016 and 2011 | CT level")

summary(bike_lane_changes$Bike_lane_diff.2011.2016ct)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1377.7     0.0     0.0   262.1   189.7 14463.8

ggplot(bike_lane_changes) +
  geom_histogram(aes(as.numeric(Bike_lane_diff.by.street.2011.2016ct))) + 
  xlab("Difference of bike lane length between 2016 and 2011, normalized by street")

summary(bike_lane_changes$Bike_lane_diff.by.street.2011.2016ct)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -100.000    0.000    0.000    3.665    4.496  100.000      252

2.2 Canopy changes

Canopy changes is based on data produced by CMM, using multispectral aerial imagery and lidar. In order to sync the observations with the census years, we focus on 2011 and 2019 with one extra observation point in 2015.

The processing steps are similar to the ones for the bike lanes:

• Import the raster for each of the 3 years
• Compute proportion of canopy within each CT for the 3 years

# Codes du raster "espace vert"
# 0. No data (hors CMM)
# 1. NDVI < 0,3 et MNH < 3,0m = Minéral bas (route, stationnement, etc.)
# 2. NDVI < 0,3 et MNH ≥ 3,0m = Minéral haut (constructions)
# 3. NDVI ≥ 0,3 et MNH < 3,0m = Végétal bas (culture, gazon, etc.)
# 4. NDVI ≥ 0,3 et MNH ≥ 3,0m = Végétal haut (canopée)
# 5. Aquatique

# Load rasters into pg database for further processing
system('psql -d gentrif_bei -c "CREATE EXTENSION IF NOT EXISTS postgis"')
system('psql -d gentrif_bei -c "CREATE EXTENSION IF NOT EXISTS postgis_raster"')

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2019') IS NOT NULL;")) == 0) {
  system("raster2pgsql -s 32188 -I -C -M data/canopy/2019/*.tif -F -t 1000x1000 canopee2019 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2019' already imported") }

## PG Raster 'canopee2019' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2017') IS NOT NULL;")) == 0) {
  system("raster2pgsql -s 32188 -I -C -M data/canopy/2017/*.tif -F -t 1000x1000 canopee2017 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2017' already imported") }

## PG Raster 'canopee2017' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2015') IS NOT NULL;")) == 0) {
  system("raster2pgsql -s 32188 -I -C -M data/canopy/2015/*.tif -F -t 1000x1000 canopee2015 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2015' already imported") }

## PG Raster 'canopee2015' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2011') IS NOT NULL;")) == 0) {
  system("raster2pgsql -s 32188 -I -C -M data/canopy/2011/*.tif -F -t 1000x1000 canopee2011 | psql -d gentrif_bei", intern = TRUE)
} else { message("PG Raster 'canopee2011' already imported") }

## PG Raster 'canopee2011' already imported

# Resample to 10m as the original rasters have a 1m resolution, which is too high to allow for a swift processing
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2019_10m') IS NOT NULL;")) == 0) {
  system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2019 mode=2\" -r mode -tr 10 10 data/canopy/canopee2019_10m.tif")
  system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2019_10m.tif -F -t 100x100 canopee2019_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2019_10m' already imported") }

## PG Raster 'canopee2019_10m' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2017_10m') IS NOT NULL;")) == 0) {
  system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2017 mode=2\" -r mode -tr 10 10 data/canopy/canopee2017_10m.tif")
  system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2017_10m.tif -F -t 100x100 canopee2017_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2017_10m' already imported") }

## PG Raster 'canopee2017_10m' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2015_10m') IS NOT NULL;")) == 0) {
  system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2015 mode=2\" -r mode -tr 10 10 data/canopy/canopee2015_10m.tif")
  system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2015_10m.tif -F -t 100x100 canopee2015_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2015_10m' already imported") }

## PG Raster 'canopee2015_10m' already imported

if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('canopee2011_10m') IS NOT NULL;")) == 0) {
  system("gdal_translate -of GTiff PG:\"host=localhost dbname=gentrif_bei table=canopee2011 mode=2\" -r mode -tr 10 10 data/canopy/canopee2011_10m.tif")
  system("raster2pgsql -s 32188 -I -C -M data/canopy/canopee2011_10m.tif -F -t 100x100 canopee2011_10m | psql -d gentrif_bei")
} else { message("PG Raster 'canopee2011_10m' already imported") }

## PG Raster 'canopee2011_10m' already imported

# Push CT16 to pg
if (nrow(dbGetQuery(con_bei, "SELECT 1 test WHERE to_regclass('ct16') IS NOT NULL;")) == 0) {
  CT16 %>%
    st_transform(crs = 32188) %>%
    st_write(con_bei, "ct16",
             layer_options = c("OVERWRITE=yes", "LAUNDER=true", "SPATIAL_INDEX=gist", "GEOMETRY_NAME=geom"))
  system("psql -d gentrif_bei -c 'CREATE INDEX ON  ct16 USING gist (geometry)'")
} else { message("PG Layer CT16 already imported") }

## PG Layer CT16 already imported

WITH cnt19 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2019_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee19 AS (
    SELECT "GeoUID"
        ,round(.1*.1 * sum(total) FILTER (WHERE value in (3, 4))) AS area_esp_vert_2019 -- area expressed in hectares
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 1) AS pct_esp_vert_2019
        ,round(.1*.1 * sum(total) FILTER (WHERE value = 4)) AS area_esp_vert_high_2019
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 1) AS pct_esp_vert_high_2019
    FROM cnt19
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt17 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2017_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee17 AS (
    SELECT "GeoUID"
        ,round(.1*.1 * sum(total) FILTER (WHERE value in (3, 4))) AS area_esp_vert_2017
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 1) AS pct_esp_vert_2017
        ,round(.1*.1 * sum(total) FILTER (WHERE value = 4)) AS area_esp_vert_high_2017
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 1) AS pct_esp_vert_high_2017
    FROM cnt17
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt15 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2015_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee15 AS (
    SELECT "GeoUID"
        ,round(.1*.1 * sum(total) FILTER (WHERE value in (3, 4))) AS area_esp_vert_2015
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 1) AS pct_esp_vert_2015
        ,round(.1*.1 * sum(total) FILTER (WHERE value = 4)) AS area_esp_vert_high_2015
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 1) AS pct_esp_vert_high_2015
    FROM cnt15
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
),
cnt11 AS (
    SELECT "GeoUID", "Population"
        ,(pvc).value, SUM((pvc).count) As total
    FROM (SELECT "GeoUID", "Population"
            ,ST_ValueCount(ST_Clip(rast, geometry)) As pvc
        FROM canopee2011_10m
        JOIN ct16 ON ST_Intersects(geometry, rast)
    ) As foo
    GROUP BY "GeoUID", "Population", (pvc).value
),
canopee11 AS (
    SELECT "GeoUID"
        ,round(.1*.1 * sum(total) FILTER (WHERE value in (3, 4))) AS area_esp_vert_2011
        ,round(100. * sum(total) FILTER (WHERE value in (3, 4)) / sum(total), 1) AS pct_esp_vert_2011
        ,round(.1*.1 * sum(total) FILTER (WHERE value = 4)) AS area_esp_vert_high_2011
        ,round(100. * sum(total) FILTER (WHERE value = 4) / sum(total), 1) AS pct_esp_vert_high_2011
    FROM cnt11
    WHERE value > 0 -- discard no data, including postgis raster no data
    GROUP BY "GeoUID", "Population"
)
SELECT "GeoUID"
    ,round((st_area(geometry) / 10000)::numeric, 1) ct_area_h -- CT area in hectares
    ,COALESCE(area_esp_vert_2011, 0) area_esp_vert_2011
    ,coalesce(area_esp_vert_high_2011, 0) area_esp_vert_high_2011
    ,coalesce(pct_esp_vert_high_2011, 0) pct_esp_vert_high_2011
    ,coalesce(pct_esp_vert_2011, 0) pct_esp_vert_2011
    ,COALESCE(area_esp_vert_2015, 0) area_esp_vert_2015
    ,coalesce(area_esp_vert_high_2015, 0) area_esp_vert_high_2015
    ,coalesce(pct_esp_vert_high_2015, 0) pct_esp_vert_high_2015
    ,coalesce(pct_esp_vert_2015, 0) pct_esp_vert_2015
    ,COALESCE(area_esp_vert_2017, 0) area_esp_vert_2017
    ,coalesce(area_esp_vert_high_2017, 0) area_esp_vert_high_2017
    ,coalesce(pct_esp_vert_high_2017, 0) pct_esp_vert_high_2017
    ,coalesce(pct_esp_vert_2017, 0) pct_esp_vert_2017
    ,COALESCE(area_esp_vert_2019, 0) area_esp_vert_2019
    ,coalesce(area_esp_vert_high_2019, 0) area_esp_vert_high_2019
    ,coalesce(pct_esp_vert_high_2019, 0) pct_esp_vert_high_2019
    ,coalesce(pct_esp_vert_2019, 0) pct_esp_vert_2019
FROM ct16
FULL JOIN canopee19 USING ("GeoUID")
FULL JOIN canopee17 USING ("GeoUID")
FULL JOIN canopee15 USING ("GeoUID")
FULL JOIN canopee11 USING ("GeoUID");

2.3 Pampalon index

Get it here

pampalon <- read.xlsx("data/Canada2016Pampalon/A-MSDIData_Can2016_eng/1. EquivalenceTableCanada2016_ENG.xlsx", sheet = 2) %>%
  mutate(DA = as.character(DA)) %>%
  select(DA, SCOREMAT, SCORESOC)

# 2016 DA boundaries for Montreal
DA16 <- get_census(dataset='CA16', regions=list(CMA='24462'), level='DA', geo_format = "sf") %>%
  filter(Type == "DA") %>%
  st_transform(st_crs(bike_lane))

## Reading geo data from local cache.

pampalon <- DA16 %>%
  inner_join(pampalon, by = c("GeoUID" = "DA")) %>%
  as.data.frame()

# Get Pampalon 2006
pampalon06 <- read.xlsx("data/Canada2006Pampalon/A-MSDIData_Can2006_eng/1. CorrespondenceTable_Can2006_eng.xlsx", sheet = 2) %>%
  mutate(DA = as.character(DA)) %>%
  select(DA, DAPOP2006, SCOREMAT, SCORESOC)

# Get LUT DA2006 <-> DA2011 from StatCan
lut_da.1 <- read.csv("data/2011_92-156_DA_AD_txt/2011_92-156_DA_AD.txt", colClasses = "character", 
                     header = FALSE, col.names = c("DAUID2011.ADIDU2011", "DAUID2006.ADIDU2006", "DBUID2011", "DA_rel_flag")) %>%
  select(!c(DBUID2011, DA_rel_flag)) %>%
  unique()

# Link Pampalon 2011 to LUT and compute weighted mean of scores of Pampalon 2011
# NB: population numbers will diverge from  reality when more than one DA is merged into one DA of next census
pampalon06.11 <- pampalon06 %>%
  inner_join(lut_da.1, by = c("DA" = "DAUID2006.ADIDU2006")) %>%
  group_by(DAUID2011.ADIDU2011) %>%
  summarise(pop2006 = sum(DAPOP2006),
            SCOREMAT.06 = weighted.mean(SCOREMAT, DAPOP2006, na.rm = TRUE),
            SCORESOC.06 = weighted.mean(SCORESOC, DAPOP2006, na.rm = TRUE))

# Get Pampalon 2011
pampalon11 <- read.xlsx("data/Canada2011Pampalon/A-MSDIData_Can2011_eng/1. CorrespondenceTable_Can2011_eng.xlsx", sheet = 2) %>%
  mutate(DA = as.character(DA)) %>%
  select(DA, DAPOP2011, SCOREMAT, SCORESOC)

# Get LUT DA2011 <-> DA2016 from StatCan
lut_da <- read.csv("data/2016_92-156_DA_AD_csv/2016_92-156_DA_AD.csv", colClasses = "character") %>%
  select(!c(DBUID2016.IDIDU2016, DA_rel_flag.AD_ind_rel)) %>%
  unique()

# Link Pampalon 2011 to LUT, then to Pampalon 06 and finally compute weighted mean of scores of Pampalon 2011
pampalon11.16 <- pampalon11 %>%
  inner_join(lut_da, by = c("DA" = "DAUID2011.ADIDU2011")) %>%
  left_join(pampalon06.11, by =c("DA" = "DAUID2011.ADIDU2011")) %>%
  group_by(DAUID2016.ADIDU2016) %>%
  summarise(pop2011 = sum(DAPOP2011),
            SCOREMAT = weighted.mean(SCOREMAT, DAPOP2011, na.rm = TRUE),
            SCORESOC = weighted.mean(SCORESOC, DAPOP2011, na.rm = TRUE),
            SCOREMAT.06 = weighted.mean(SCOREMAT.06, pop2006, na.rm = TRUE),
            SCORESOC.06 = weighted.mean(SCORESOC.06, pop2006, na.rm = TRUE),
            pop2006 = sum(pop2006))

# Then link Pampalon 2011 to 2016
pampalon <- pampalon %>%
  left_join(pampalon11.16, by = c("GeoUID" = "DAUID2016.ADIDU2016"), suffix = c(".16", ".11"))

# Aggregate at the CT level
pampalon_CT <- pampalon %>%
  group_by(CT_UID) %>%
  summarise(wSCOREMAT.2016 = weighted.mean(SCOREMAT.16, Population, na.rm = TRUE),
            wSCORESOC.2016 = weighted.mean(SCORESOC.16, Population, na.rm = TRUE),
            wSCOREMAT.2011 = weighted.mean(SCOREMAT.11, pop2011, na.rm = TRUE),
            wSCORESOC.2011 = weighted.mean(SCORESOC.11, pop2011, na.rm = TRUE),
            wSCOREMAT.2006 = weighted.mean(SCOREMAT.06, pop2006, na.rm = TRUE),
            wSCORESOC.2006 = weighted.mean(SCORESOC.06, pop2006, na.rm = TRUE))

# Clean up
rm(lut_da, lut_da.1, pampalon11.16, pampalon06.11, pampalon11, pampalon06)

# Display map
.pampalon_CT_geom <- CT16 %>%
  left_join(pampalon_CT, by = c("GeoUID" = "CT_UID")) %>%
  filter(interact_aoi)

.pampalon_data <- bi_class(.pampalon_CT_geom, x = wSCOREMAT.2016, y = wSCORESOC.2016, style = "quantile", dim = 3)

## Warning in classInt::classIntervals(bins_x, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes

## Warning in classInt::classIntervals(bins_y, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes

.map <- ggplot() + 
  geom_sf(data = .pampalon_data, mapping = aes(fill = bi_class), color = "white", size = 0.1, show.legend = FALSE) +
  bi_scale_fill(pal = "DkBlue", dim = 3) +
  labs(title = "Pampalon: material and social deprivation index") + 
  theme(panel.background = element_rect(fill = "white"),
        #axis.ticks = element_blank(),
        #axis.text = element_blank(),
        panel.grid = element_line(color = "darkgray", size = 0.2))
.legend <- bi_legend(pal = "DkBlue",
                    dim = 3,
                    xlab = "Material ",
                    ylab = "Social ",
                    size = 8)
ggdraw() +
  draw_plot(.map, 0, 0, 1, 1) +
  draw_plot(.legend, 0.1, .7, 0.2, 0.2)

2.4 Gentrification

Using Ding metric computed on 5 year span.

# Load gentrified CTs, 5 year span (from repo gentrification_metrics)
ding <- list()
ding[["2016"]] <- st_read("data/gentrified_5years.gpkg", "gentrified_ding_16", quiet=TRUE) %>%
  filter(cma_uid_16 == "24462") %>%
  st_transform(st_crs(bike_lane))
ding[["2011"]] <- st_read("data/gentrified_5years.gpkg", "gentrified_ding_11", quiet=TRUE) %>%
  filter(cma_uid_11 == "24462") %>%
  st_transform(st_crs(bike_lane))
ding[["2006"]] <- st_read("data/gentrified_5years.gpkg", "gentrified_ding_06", quiet=TRUE) %>%
  filter(cma_uid_06 == "24462") %>%
  st_transform(st_crs(bike_lane))

.ding_map <- ding[["2016"]] %>%
  left_join(select(as.data.frame(CT16), GeoUID, interact_aoi), by = c("ct_uid_16" = "GeoUID")) %>%
  filter(interact_aoi)

ggplot(data = .ding_map) + 
  geom_sf(aes(fill = gentrified_2016_2011, colour=gentrifiable_2011)) +
  scale_fill_manual(values = c("gray", "red", "darkgray"), name = "Gentrified in 2016") +
  scale_colour_manual(values = c("darkgray", "darkred", "darkgray"), name = "Gentrifiable in 2011") +
  labs(title = "Census tract gentrification status in 2016")

2.5 Equity metrics | low-income population and visible minority

Introduced here as a proposition, nothing acted (2022-02-04)

# Visible Minority
# - v_CA16_3954: Total - Visible minority for the population in private households - 25% sample data (Total)
# - v_CA16_3957: Total visible minority population (Total)

# Low income (LIM-AT)
# - v_CA16_2540: Prevalence of low income based on the Low-income measure, after tax (LIM-AT) (%) (Total)
equity_ct16 <- get_census(dataset='CA16', regions=list(CMA='24462'), level='CT', geo_format = "sf",
                          vectors = c("v_CA16_3954", "v_CA16_3957", "v_CA16_2540")) %>%
  filter(Type == "CT") %>%
  transmute(CT_UID = GeoUID,
            vis_minority_2016 = `v_CA16_3957: Total visible minority population` / `v_CA16_3954: Total - Visible minority for the population in private households - 25% sample data` * 100,
            low_income_2016 = `v_CA16_2540: Prevalence of low income based on the Low-income measure, after tax (LIM-AT) (%)`)

## Reading vectors data from local cache.

## Reading geo data from local cache.

# Visible Minority
# - v_CA11N_457: CA 2011 NHS, Total population in private households by visible minority (Total)
# - v_CA11N_460: CA 2011 NHS, Total population in private households by visible minority, Total visible minority population (Total)

# Low income (LIM-AT)
# - v_CA11N_2606: CA 2011 NHS, Prevalence of low income in 2010 based on after-tax low-income measure % (Total)
equity_ct11 <- get_census(dataset='CA11', regions=list(CMA='24462'), level='CT', geo_format = "sf",
                          vectors = c("v_CA11N_457", "v_CA11N_460", "v_CA11N_2606")) %>%
  filter(Type == "CT") %>%
  transmute(CT_UID = GeoUID,
            vis_minority_2011 = `v_CA11N_460: Total visible minority population` / `v_CA11N_457: Total population in private households by visible minority` * 100,
            low_income_2011 = `v_CA11N_2606: Prevalence of low income in 2010 based on after-tax low-income measure %`)

## Reading vectors data from local cache.
## Reading geo data from local cache.

# Visible Minority
# - v_CA06_1302: Total population by visible minority groups
# - v_CA06_1303: Total population by visible minority groups, Total visible minority population

# Low income (LIM-AT)
# - v_TX2006_551: After-tax low income status of tax filers and dependents (census family low income measure, CFLIM-AT) for couple and lone parent families by family composition, 2006 | All family units | Persons in Low Income | % - Total
equity_ct06 <- get_census(dataset='CA06', regions=list(CMA='24462'), level='CT', geo_format = "sf",
                          vectors = c("v_CA06_1302", "v_CA06_1303", "v_TX2006_551")) %>%
  filter(Type == "CT") %>%
  transmute(CT_UID = GeoUID,
            vis_minority_2006 = `v_CA06_1303: Total visible minority population` / `v_CA06_1302: Total population by visible minority groups - 20% sample data` * 100,
            low_income_2006 = `v_TX2006_551: % - Total`)

## Reading vectors data from local cache.
## Reading geo data from local cache.

equity_ct <- st_join(equity_ct16, equity_ct11, left=TRUE, largest=TRUE, suffix=c("", "_2011")) %>% # join on largest overlap, to overcome mismatch in CT UID
  st_join(equity_ct06, left=TRUE, largest=TRUE, suffix=c("", "_2006")) %>%
  data.frame()

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

# cleanup
rm(equity_ct11, equity_ct16, equity_ct06)

# Display map
.equity_CT_geom <- CT16 %>%
  left_join(equity_ct, by = c("GeoUID" = "CT_UID")) %>%
  filter(interact_aoi)

.equity_data <- bi_class(.equity_CT_geom, x = vis_minority_2016, y = low_income_2016, style = "quantile", dim = 3)

## Warning in classInt::classIntervals(bins_x, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes

## Warning in classInt::classIntervals(bins_y, n = dim, style = "quantile"): var
## has missing values, omitted in finding classes

.map <- ggplot() + 
  geom_sf(data = .equity_data, mapping = aes(fill = bi_class), color = "white", size = 0.1, show.legend = FALSE) +
  bi_scale_fill(pal = "Brown", dim = 3) +
  labs(title = "Equity metrics: % of visible minority and % of low-income household") + 
  theme(panel.background = element_rect(fill = "white"),
        #axis.ticks = element_blank(),
        #axis.text = element_blank(),
        panel.grid = element_line(color = "darkgray", size = 0.2))
.legend <- bi_legend(pal = "Brown",
                    dim = 3,
                    xlab = "Vis. Minority ",
                    ylab = "Low-Income ",
                    size = 8)
ggdraw() +
  draw_plot(.map, 0, 0, 1, 1) +
  draw_plot(.legend, 0.1, .7, 0.2, 0.2)

2.6 Build complete dataset

All variables + outcome linked at the CT level + quintiles of SES variables

.bike_lane_changes <- bike_lane_changes %>%
  as.data.frame() %>%
  select(GeoUID, starts_with("street_length"), ends_with("ct", ignore.case = FALSE)) %>%
  select(GeoUID, starts_with("street_length"), starts_with("Bike_lane")) # Drop individual category lane length

bei_df <- CT16 %>%
  as.data.frame() %>%
  transmute(CT_UID = GeoUID,
            CD_UID = CD_UID,
            CSD_UID = CSD_UID,
            interact_aoi = interact_aoi,
            zone = case_when(CD_UID == "2466" ~ "Montreal",
                             CD_UID == "2458" ~ "Longueuil",
                             CD_UID == "2465" ~ "Laval",
                             TRUE ~ "Other"),
            Population = Population) %>%
  left_join(pampalon_CT, by="CT_UID") %>%
  left_join(select(as.data.frame(ding$`2016`), ct_uid_16, starts_with("gentrif")), by=c("CT_UID" = "ct_uid_16")) %>%
  left_join(select(as.data.frame(ding$`2011`), ct_uid_11, starts_with("gentrif")), by=c("CT_UID" = "ct_uid_11")) %>%
  left_join(select(as.data.frame(ding$`2006`), ct_uid_06, starts_with("gentrif")), by=c("CT_UID" = "ct_uid_06")) %>%
  left_join(select(as.data.frame(equity_ct), !c("geometry", "CT_UID_2011", "CT_UID_2006")), by="CT_UID") %>%
  left_join(.bike_lane_changes, by=c("CT_UID" = "GeoUID")) %>%
  left_join(as.data.frame(rename_with(esp_vert_ct, ~ paste0(., "ct"))), by=c("CT_UID" = "GeoUIDct")) %>%
  # Compute quintile of SES variables
  mutate(wSCOREMAT.2006.Q = ntile(wSCOREMAT.2006, 5),
         wSCOREMAT.2011.Q = ntile(wSCOREMAT.2011, 5),
         wSCOREMAT.2016.Q = ntile(wSCOREMAT.2016, 5),
         low_income_2006.Q = ntile(low_income_2006, 5),
         low_income_2011.Q = ntile(low_income_2011, 5),
         low_income_2016.Q = ntile(low_income_2016, 5),
         vis_minority_2006.Q = ntile(vis_minority_2006, 5),
         vis_minority_2011.Q = ntile(vis_minority_2011, 5),
         vis_minority_2016.Q = ntile(vis_minority_2016, 5)) %>%
  units::drop_units()
  
  
head(bei_df)

write.csv(bei_df, "data/_results/bei_equity.csv", na="", row.names = FALSE)

# keep only interact CT
bei_df_aoi <- filter(bei_df, interact_aoi)

Included variables:

• Census Tracts variables
  • CT_UID: 2016 Census Tract ID
  • CD_UID: 2016 Census Division
  • CSD_UID: 2016 Census Subdivision
  • interact_aoi: Does CT belong to INTERACT study area?
  • zone: Montreal / Laval / Longueuil, etc. / Other
  • Population: 2016 Population within CT
  • ct_area_m2ct: Area of CT, in square meters
• Gentrification metrics
  • gentrified_2016_2011: Is the CT gentrified in 2016?
  • gentrifiable_2011: Is the CT candidate to gentrification in 2011?
  • gentrified_2011_2006: Is the CT gentrified in 2011
  • gentrifiable_2006: Is the CT candidate to gentrification in 2006
  • gentrified_2006_2001: Is the CT gentrified in 2006
  • gentrifiable_2001: Is the CT candidate to gentrification in 2001
• Pampalon’s metrics
  • wSCORESOC.2016: Social deprivation index in 2016 (population weighted)
  • wSCOREMAT.2016: Material deprivation index in 2016 (population weighted)
  • wSCOREMAT.2016.Q: Quintile of material deprivation index in 2016 (population weighted)
  • wSCORESOC.2011: Social deprivation index in 2011 (population weighted)
  • wSCOREMAT.2011: Material deprivation index in 2011 (population weighted)
  • wSCOREMAT.2011.Q: Quintile of aterial deprivation index in 2011 (population weighted)
  • wSCORESOC.2006: Social deprivation index in 2006 (population weighted)
  • wSCOREMAT.2006: Material deprivation index in 2006 (population weighted)
  • wSCOREMAT.2006.Q: Quintile of material deprivation index in 2006 (population weighted)
• Social profile
  • vis_minority_2016: % of visible minority in CT 2016
  • vis_minority_2016.Q: Quintile of % of visible minority in CT 2016
  • low_income_2016: prevalence of low income in CT 2016
  • low_income_2016.Q: Quintile of prevalence of low income in CT 2016
  • vis_minority_2011: % of visible minority in CT 2011
  • vis_minority_2011.Q: Quintile of % of visible minority in CT 2011
  • low_income_2011: prevalence of low income in CT 2011
  • low_income_2011.Q: Quintile of prevalence of low income in CT 2011
  • vis_minority_2006: % of visible minority in CT 2006
  • vis_minority_2006.Q: Quintile of % of visible minority in CT 2006
  • low_income_2006: prevalence of low income in CT 2006
  • low_income_2006.Q: Qintile of prevalence of low income in CT 2006
• Bike lane density
  • Bike_lane_total.{2016|2011}ct: total length of bike lanes, in 2016 or 2011, within CT
  • Bike_lane_total.hm.{2016|2011}ct: total length of bike lanes, in hectometers (as integers), in 2016 or 2011, within CT
  • Bike_lane.by.street.{2016|2011}ct: % of bike lanes compared to streets, in 2016 or 2011, within CT
  • Bike_lane_diff.2011.2016ct: change in total length of bike lanes between 2011 and 2016, within CT
  • Bike_lane_diff.by.street.2011.2016ct: change in total length of bike lanes between 2011 and 2011, normalized by street length, within CT
  • Bike_lane_diff.by.area.2011.2016ct: change in total length of bike lanes between 2011 and 2011, normalized by area, within CT
• Green spaces
  • area_esp_vert_{2011|2015|2019}ct: area of green space in 2011, 2015 or 2019 within CT (in hectares)
  • area_esp_vert_high_{2011|2015|2019}ct: same as above, except for trees (high canopy)
  • pct_esp_vert_{2011|2015|2019}ct: % of green space in 2011, 2015 or 2019 within CT
  • pct_esp_vert_high_{2011|2015|2019}ct: same as above, except for trees (high canopy)
  • pct_esp_vert_diff{2011|2015}.{2015|2019}ct: change in % of green space between 2011 and 2015, 2011 and 2019 as well as 2011 and 2019, within CT

3 Preliminary analyses

INTERACT study area ~ Montréal, Laval, Longueuil, Brossard, St-Lambert

3.1 SES variable distribution

.bei_df_long <- bei_df_aoi %>% 
  select(CT_UID, CD_UID, starts_with("wSCORE")) %>%
  select(!ends_with('.Q')) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_histogram() + 
  facet_wrap(~name) #, scales = "free")

## Warning: Removed 86 rows containing non-finite values (stat_bin).

.bei_df_long <- bei_df_aoi %>%
  select(CT_UID, CD_UID, starts_with("vis_minority"), starts_with("low_income")) %>%
  select(!ends_with('.Q')) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_histogram() + 
  facet_wrap(~name) #, scales = "free")

## Warning: Removed 85 rows containing non-finite values (stat_bin).

.bei_df_long <- bei_df_aoi %>%
  select(CT_UID, CD_UID, starts_with("gentrif")) %>%
  select(!ends_with('.Q')) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_bar() + 
  facet_wrap(~name) #, scales = "free", ncol = 3)

3.2 BEI variable distributions

3.2.1 Absolute values

.bei_df_long <- bei_df_aoi %>% 
  filter(interact_aoi) %>%
  select(CT_UID, CD_UID, matches("^Bike_lane_total.*ct$"), matches("^area_esp_vert.*ct$")) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_histogram() + 
  facet_wrap(~name, scales = "free", ncol = 4)

3.2.2 Relative values

.bei_df_long <- bei_df_aoi %>% 
  filter(interact_aoi) %>%
  select(CT_UID, CD_UID, matches("^Bike_lane.by.*ct$"), matches("^pct_esp_vert.*ct$")) %>%
  pivot_longer(!c(CT_UID, CD_UID))

ggplot(.bei_df_long, aes(value)) +
  geom_histogram() + 
  facet_wrap(~name, scales = "free", ncol = 3)

4 Association between SES and Urban Conditions at baseline (2011)

Looking at objective #1: do urban interventions tend to be located in low SES neighborhoods?. We look at \[Urban Condition_{2011} = f(SES_{2011})\] as well as \[Urban Condition_{2011} = f(Gentrification_{2011 \to 2016})\]

Here \(UrbanCondition\) means the state of the urban environment features, such as length of bike lanes, greenness coverage, etc. at one specific moment. This needs to be distinguished from \(UrbanIntervention\), which accounts for the changes in the \(UrbanConditions\) between two years (see below).

We fit LMM models, with spatial random effect using spaMM package.

4.1 UC vs Pampalon | material

4.1.1 Bike lane length

Bike lane ratio to streets at the Census tract level. We use an offset, as explained in Anderson’s presentation, to handle the denominator and model the Bike lane length in hectometers using a Poisson distribution.

f <- Bike_lane_total.hm.2011ct ~ wSCOREMAT.2011.Q + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(Bike_lane_total.hm.2011ct ~ wSCOREMAT.2011.Q + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2011ct ~ wSCOREMAT.2011.Q + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                  Estimate Cond. SE t-value   p-value
## (Intercept)        -2.387  0.18193 -13.118 0.0000000
## wSCOREMAT.2011.Q   -0.154  0.04662  -3.303 0.0009572
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1369437 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  1.532  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2262.435

ci.lmm <- confint(res.lmm, 'wSCOREMAT.2011.Q')

## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q 
##            -0.24668125            -0.06222581

# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - wSCOREMAT.2011.Q)
LRT(res.lmm.0, res.lmm)

##      chi2_LR df     p_value
## p_v 10.77208  1 0.001030424

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# AIC
AIC(res.lmm)

Interpretation of the wSCOREMAT.2011.Q coefficient (significant): we have \(e^{-0.154}\) = 0.857 (95%CI: 0.781 – 0.94), which means in CTs with bike lanes, each increase of 1 quintile in the Material Score 2011 leads to a 14.3% decrease the ratio of bike lanes to street in 2011.

4.1.2 Bike lane length (hurdle)

Resorting to a hurdle model, to accommodate the zero-inflated distribution.

f <- Bike_lane_total.hm.2011ct ~ wSCOREMAT.2011.Q + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

# Zero model
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2011ct != 0) ~ wSCOREMAT.2011.Q + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))

## formula: I(Bike_lane_total.hm.2011ct != 0) ~ wSCOREMAT.2011.Q + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                  Estimate Cond. SE t-value   p-value
## (Intercept)       -2.4427  0.35400  -6.900 5.186e-12
## wSCOREMAT.2011.Q  -0.2107  0.08799  -2.395 1.662e-02
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1377454 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  1.461  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -354.6466

ci.lmm.0 <- confint(res.lmm.zero, 'wSCOREMAT.2011.Q')

## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q 
##            -0.40742595            -0.03286067

res.lmm.0 <- update(res.lmm.zero, . ~ . -wSCOREMAT.2011.Q)
LRT(res.lmm.0, res.lmm.zero)

##     chi2_LR df    p_value
## p_v 5.41559  1 0.01995771

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.zero, plot = F)
  plot(simout, title = "DHARMa residual | Logit")
}

# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
  filter(Bike_lane_total.hm.2011ct != 0) %>%
  tidy_df(CT16, f)

res.lmm.cnt <- fitme(Bike_lane_total.hm.2011ct ~ wSCOREMAT.2011.Q + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2011ct ~ wSCOREMAT.2011.Q + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                  Estimate Cond. SE t-value p-value
## (Intercept)      -2.03822  0.09322 -21.864 0.00000
## wSCOREMAT.2011.Q -0.06074  0.02640  -2.301 0.02139
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1338634 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.3906  
## # of obs: 469; # of groups: ct_no, 469 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1735.449

ci.lmm.cnt <- confint(res.lmm.cnt, 'wSCOREMAT.2011.Q')

## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q 
##           -0.113107222           -0.008677374

res.lmm.0 <- update(res.lmm.cnt, . ~ . -wSCOREMAT.2011.Q)
LRT(res.lmm.0, res.lmm.cnt)

##      chi2_LR df    p_value
## p_v 5.223183  1 0.02228772

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout, title = "DHARMa residual | Count (truncated Poisson)")
}

# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)

# Combined AIC
"Combined Hurdle AIC:"

## [1] "Combined Hurdle AIC:"

-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)

## [1] 4188.19

Interpretation of the wSCOREMAT.2011.Q coefficient (significant in both sub-models): - we have \(e^{-0.211}\) = 0.81 (95%CI: 0.665 – 0.968), which means each increase of 1 quintile in the Material Score 2011 leads to a 19% decrease the chance of having a bicycle lane in the CT - we have \(e^{-0.061}\) = 0.941 (95%CI: 0.893 – 0.991), which means in CTs with bike lanes, each increase of 1 quintile in the Material Score 2011 leads to a 5.9% decrease the ratio of bike lanes to street in 2011.

4.1.3 Greenness

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 at the Census tract level, using a Poisson distribution

f <- area_esp_vert_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx, 
                 family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                  Estimate Cond. SE t-value   p-value
## (Intercept)      -0.62343  0.03922 -15.895 0.000e+00
## wSCOREMAT.2011.Q -0.08301  0.01029  -8.065 7.772e-16
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1383759 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.0557  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2413.744

ci.lmm <- confint(res.lmm, 'wSCOREMAT.2011.Q')

## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q 
##            -0.10320178            -0.06275188

# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - wSCOREMAT.2011.Q)
LRT(res.lmm.0, res.lmm)

##      chi2_LR df      p_value
## p_v 60.94647  1 5.884182e-15

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

Interpretation of the wSCOREMAT.2011.Q coefficient (significant): we have \(e^{-0.083}\) = 0.92 (95%CI: 0.902 – 0.939), which means each increase of 1 quintile in the Material Score leads to a 8% decrease in the percentage of greenspace area within CT in 2011.

4.1.4 Greenness (Neg Bin)

No better than Poisson

f <- area_esp_vert_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

# Negbin
res.lmm.negbin <- fitme(area_esp_vert_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = negbin())
summary(res.lmm.negbin, details = c(p_value="Wald"))

## formula: area_esp_vert_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars, lambda and NB_shape by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda and NB_shape by 'outer' ML, maximizing p_v.
## family: Neg.binomial(shape=816000)( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                  Estimate Cond. SE t-value   p-value
## (Intercept)      -0.62344  0.03922 -15.895 0.000e+00
## wSCOREMAT.2011.Q -0.08301  0.01029  -8.065 7.772e-16
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1383752 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.0557  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2413.745

#ci.lmm <- confint(res.lmm.negbin, 'wSCOREMAT.2011.Q')

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.negbin, plot = F)
  plot(simout, title = "DHARMa residual | Count, negative binomial")
}

4.1.5 Greenness (hurdle)

Not run due to complete separation error in zero model:

The following terms are causing separation among the sample points: (Intercept), wSCOREMAT.2011.Q 
    some estimates of fixed-effect coefficients could be practically infinite, 
    causing numerical issues in various functions.

# f <- area_esp_vert_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct))
# 
# clean_bei <- tidy_df(bei_df_aoi, CT16, f)

# # Zero model
# res.lmm.zero <- fitme(I(area_esp_vert_2011ct != 0) ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) + adjacency(1|ct_no),
#                  data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
#                  family = binomial())
# summary(res.lmm.zero, details = c(p_value="Wald"))
# 
# ci.lmm.0 <- confint(res.lmm.zero, 'wSCOREMAT.2011.Q')
# 
# res.lmm.0 <- update(res.lmm.zero, . ~ . -wSCOREMAT.2011.Q)
# LRT(res.lmm.0, res.lmm.zero)

# if (plot.resid) {
#   simout <- simulateResiduals(fittedModel = res.lmm.zero, plot = F)
#   plot(simout, title = "DHARMa residual | Logit")
# }
# 
# # Count model (Poisson)
# clean_bei.cnt <- bei_df_aoi %>%
#   filter(area_esp_vert_2011ct != 0) %>%
#   tidy_df(CT16, f)
# 
# res.lmm.cnt <- fitme(area_esp_vert_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) + adjacency(1|ct_no),
#                  data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
#                  family = Tpoisson())
# summary(res.lmm.cnt, details = c(p_value="Wald"))
# 
# ci.lmm.cnt <- confint(res.lmm.cnt, 'wSCOREMAT.2011.Q')
# 
# res.lmm.0 <- update(res.lmm.cnt, . ~ . -wSCOREMAT.2011.Q)
# LRT(res.lmm.0, res.lmm.cnt)
# 
# if (plot.resid) {
#   simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
#   plot(simout, title = "DHARMa residual | Count (truncated Poisson)")
# }

4.1.6 Canopy (trees)

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 at the Census tract level, using a Poisson distribution

f <- area_esp_vert_high_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_high_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx, 
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_high_2011ct ~ wSCOREMAT.2011.Q + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                  Estimate Cond. SE t-value   p-value
## (Intercept)       -1.3028  0.05880 -22.155 0.000e+00
## wSCOREMAT.2011.Q  -0.1204  0.01551  -7.764 8.216e-15
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1376817 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.1304  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2154.547

ci.lmm <- confint(res.lmm, 'wSCOREMAT.2011.Q')

## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q 
##            -0.15088539            -0.08990383

# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - wSCOREMAT.2011.Q)
LRT(res.lmm.0, res.lmm)

##     chi2_LR df      p_value
## p_v 57.1245  1 4.085621e-14

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

Interpretation of the wSCOREMAT.2011.Q coefficient (significant): we have \(e^{-0.12}\) = 0.887 (95%CI: 0.86 – 0.914), which means each increase of 1 quintile in the Material Score leads to a 11.3% decrease in the percentage of tree canopy area within CT in 2011.

4.2 UC vs Visible minority

4.2.1 Bike lane length

Bike lane ratio to streets (in %) at the Census tract level

f <- Bike_lane_total.hm.2011ct ~ vis_minority_2011.Q + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(Bike_lane_total.hm.2011ct ~ vis_minority_2011.Q + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2011ct ~ vis_minority_2011.Q + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                     Estimate Cond. SE t-value  p-value
## (Intercept)          -2.3555  0.23166 -10.168 0.000000
## vis_minority_2011.Q  -0.1555  0.06033  -2.578 0.009952
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1363885 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  1.539  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2261.415

ci.lmm <- confint(res.lmm, 'vis_minority_2011.Q')

## lower vis_minority_2011.Q upper vis_minority_2011.Q 
##                -0.2755232                -0.0358404

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# AIC
AIC(res.lmm)

Interpretation of the wSCOREMAT.2011.Q coefficient: we have \(e^{NA}\) = NA (95%CI: 0.759 – 0.965), which means in CTs with bike lanes, each increase of 1 quintile in the Material Score 2011 leads to a NA% decrease the ratio of bike lanes to street in 2011.

4.2.2 Bike lane length (hurdle)

Bike lane ratio to streets (in %) at the Census tract level

f <- Bike_lane_total.hm.2011ct ~ vis_minority_2011.Q + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

# Zero model
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2011ct != 0) ~ vis_minority_2011.Q + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))

## formula: I(Bike_lane_total.hm.2011ct != 0) ~ vis_minority_2011.Q + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                     Estimate Cond. SE t-value   p-value
## (Intercept)           -2.394   0.4424  -5.411 6.255e-08
## vis_minority_2011.Q   -0.221   0.1107  -1.995 4.600e-02
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1377773 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  1.351  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -355.6222

ci.lmm.0 <- confint(res.lmm.zero, 'vis_minority_2011.Q')

## lower vis_minority_2011.Q upper vis_minority_2011.Q 
##              -0.453408008               0.003914703

res.lmm.0 <- update(res.lmm.zero, . ~ . - vis_minority_2011.Q)
LRT(res.lmm.0, res.lmm.zero)

##      chi2_LR df    p_value
## p_v 3.710964  1 0.05405616

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.zero, plot = F)
  plot(simout, title = "DHARMa residual | Logit")
}

# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
  filter(Bike_lane_total.hm.2011ct != 0) %>%
  tidy_df(CT16, f)

res.lmm.cnt <- fitme(Bike_lane_total.hm.2011ct ~ vis_minority_2011.Q + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2011ct ~ vis_minority_2011.Q + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                     Estimate Cond. SE t-value  p-value
## (Intercept)          -1.8531  0.12541 -14.776 0.000000
## vis_minority_2011.Q  -0.1086  0.03356  -3.236 0.001212
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1306523 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.3814  
## # of obs: 469; # of groups: ct_no, 469 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1728.153

ci.lmm.cnt <- confint(res.lmm.cnt, 'vis_minority_2011.Q')

## lower vis_minority_2011.Q upper vis_minority_2011.Q 
##               -0.17498489               -0.04251278

res.lmm.0 <- update(res.lmm.cnt, . ~ . - vis_minority_2011.Q)
LRT(res.lmm.0, res.lmm.cnt)

##      chi2_LR df     p_value
## p_v 10.27256  1 0.001350232

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout, title = "DHARMa residual | Count (truncated Poisson)")
}

# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)

# Combined AIC
"Combined Hurdle AIC:"

## [1] "Combined Hurdle AIC:"

-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)

## [1] 4175.551

Interpretation of the wSCOREMAT.2011.Q coefficient: - model zero-hurdle (non significant): we have \(e^{NA}\) = NA (95%CI: 0.635 – 1.004), which means each increase of 1 quintile in the proportion of visible minority 2011 leads to a NA% decrease the chance of having a bicycle lane in the CT - model count: we have \(e^{NA}\) = NA (95%CI: 0.839 – 0.958), which means in CTs with bike lanes, each increase of 1 quintile in the Material Score 2011 leads to a NA% decrease the ratio of bike lanes to street in 2011.

4.2.3 Greenness

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx, 
                 family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                     Estimate Cond. SE t-value   p-value
## (Intercept)         -0.57028  0.05165  -11.04 0.000e+00
## vis_minority_2011.Q -0.09044  0.01356   -6.67 2.559e-11
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1385817 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.05812  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                        logLik
## p_v(h) (marginal L): -2417.47

ci.lmm <- confint(res.lmm, 'vis_minority_2011.Q')

## lower vis_minority_2011.Q upper vis_minority_2011.Q 
##               -0.11703521               -0.06375021

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

Interpretation of the vis_minority_2011.Q coefficient (significant): we have \(e^{-0.09}\) = 0.914 (95%CI: 0.89 – 0.938), which means each increase of 1 quintile in the proportion of visible minority leads to a 8.6% decrease in the percentage of greenspace area within CT in 2011.

4.2.4 Canopy (trees)

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_high_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_high_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx, 
                 family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_high_2011ct ~ vis_minority_2011.Q + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                     Estimate Cond. SE t-value   p-value
## (Intercept)          -1.2040  0.07727 -15.581 0.000e+00
## vis_minority_2011.Q  -0.1388  0.02035  -6.822 8.965e-12
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1379923 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.135  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                        logLik
## p_v(h) (marginal L): -2156.85

ci.lmm <- confint(res.lmm, 'vis_minority_2011.Q')

## lower vis_minority_2011.Q upper vis_minority_2011.Q 
##               -0.17879779               -0.09883053

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

Interpretation of the vis_minority_2011.Q coefficient (significant): we have \(e^{-0.139}\) = 0.87 (95%CI: 0.836 – 0.906), which means each increase of 1 quintile in the proportion of visible minority leads to a 13% decrease in the percentage of tree canopy area within CT in 2011.

4.3 UC vs gentrified CT

Gentrified CT between 2011 and 2016

4.3.1 Bike lane length

Bike lane ratio to streets (in %) at the Census tract level

f <- Bike_lane_total.hm.2011ct ~ gentrified_2016_2011 + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(Bike_lane_total.hm.2011ct ~ gentrified_2016_2011 + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2011ct ~ gentrified_2016_2011 + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                          Estimate Cond. SE  t-value p-value
## (Intercept)               -2.8594   0.1181 -24.2215  0.0000
## gentrified_2016_2011TRUE  -0.0663   0.1341  -0.4942  0.6211
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##    1.rho 
## 0.137036 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  1.561  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2264.528

ci.lmm <- confint(res.lmm, 'gentrified_2016_2011TRUE')

## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE 
##                     -0.3324352                      0.1985979

# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - gentrified_2016_2011)
LRT(res.lmm.0, res.lmm)

##       chi2_LR df   p_value
## p_v 0.2411384  1 0.6233851

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# AIC
AIC(res.lmm)

Interpretation of the wSCOREMAT.2011.Q coefficient: none significant - model count: we have \(e^{NA}\) = NA (95%CI: 0.717 – 1.22)

4.3.2 Bike lane length (hurdle)

Bike lane ratio to streets (in %) at the Census tract level

f <- Bike_lane_total.hm.2011ct ~ gentrified_2016_2011 + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

# Zero model
res.lmm.zero <- fitme(I(Bike_lane_total.hm.2011ct != 0) ~ gentrified_2016_2011 + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))

## formula: I(Bike_lane_total.hm.2011ct != 0) ~ gentrified_2016_2011 + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                          Estimate Cond. SE t-value p-value
## (Intercept)              -3.12397   0.1952 -16.003  0.0000
## gentrified_2016_2011TRUE -0.06627   0.2473  -0.268  0.7887
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1378747 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  1.496  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -357.4465

ci.lmm.0 <- confint(res.lmm.zero, 'gentrified_2016_2011TRUE')

## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE 
##                     -0.6062718                      0.4522833

res.lmm.0 <- update(res.lmm.zero, . ~ . - gentrified_2016_2011)
LRT(res.lmm.0, res.lmm.zero)

##        chi2_LR df   p_value
## p_v 0.06253228  1 0.8025374

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.zero, plot = F)
  plot(simout, title = "DHARMa residual | Logit")
}

# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
  filter(Bike_lane_total.hm.2011ct != 0) %>%
  tidy_df(CT16, f)

res.lmm.cnt <- fitme(Bike_lane_total.hm.2011ct ~ gentrified_2016_2011 + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2011ct ~ gentrified_2016_2011 + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                          Estimate Cond. SE t-value p-value
## (Intercept)              -2.26118  0.05223 -43.290  0.0000
## gentrified_2016_2011TRUE  0.09494  0.08002   1.186  0.2355
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1267498 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.3946  
## # of obs: 469; # of groups: ct_no, 469 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1732.647

ci.lmm.cnt <- confint(res.lmm.cnt, 'gentrified_2016_2011TRUE')

## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE 
##                    -0.06921123                     0.25966646

res.lmm.0 <- update(res.lmm.cnt, . ~ . - gentrified_2016_2011)
LRT(res.lmm.0, res.lmm.cnt)

##      chi2_LR df   p_value
## p_v 1.285916  1 0.2568019

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout, title = "DHARMa residual | Count (truncated Poisson)")
}

# AIC
AIC(res.lmm.zero)
AIC(res.lmm.cnt)

# Combined AIC
"Combined Hurdle AIC:"

## [1] "Combined Hurdle AIC:"

-2 * res.lmm.zero$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)

## [1] 4188.186

Interpretation of the wSCOREMAT.2011.Q coefficient: none significant - model zero-hurdle (non significant): we have \(e^{NA}\) = NA (95%CI: 0.545 – 1.572) - model count: we have \(e^{NA}\) = NA (95%CI: 0.933 – 1.296)

4.3.3 Greenness

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx, 
                 family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                          Estimate Cond. SE t-value   p-value
## (Intercept)               -0.8304  0.02578 -32.205 0.000e+00
## gentrified_2016_2011TRUE  -0.1964  0.03343  -5.875 4.234e-09
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1385641 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.05997  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2421.862

ci.lmm <- confint(res.lmm, 'gentrified_2016_2011TRUE')

## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE 
##                     -0.2620145                     -0.1306189

# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - gentrified_2016_2011)
LRT(res.lmm.0, res.lmm)

##      chi2_LR df      p_value
## p_v 33.56416  1 6.895166e-09

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

Interpretation of the gentrified_2016_2011 coefficient (significant): we have \(e^{-0.196}\) = 0.822 (95%CI: 0.769 – 0.878), which means gentrified CTs exhibit a 17.8% decrease in the percentage of greenspace area within CT in 2011.

4.3.4 Greenness (hurdle)

Not working (complete separation)

f <- area_esp_vert_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

# Zero model
res.lmm.zero <- fitme(I(area_esp_vert_2011ct != 0) ~ gentrified_2016_2011 + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = binomial())
summary(res.lmm.zero, details = c(p_value="Wald"))

## formula: I(area_esp_vert_2011ct != 0) ~ gentrified_2016_2011 + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                          Estimate Cond. SE   t-value p-value
## (Intercept)                998.16  3063087 3.259e-04  0.9997
## gentrified_2016_2011TRUE    25.03  5570855 4.493e-06  1.0000
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##       1.rho 
## -0.06104512 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.4838  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                           logLik
## p_v(h) (marginal L): 6.82121e-13
## linear predictor estimation did not converge; try control.HLfit=list(LevenbergM=TRUE), or increase 'control.HLfit$max.iter.mean' above 200

#ci.lmm.0 <- confint(res.lmm.zero, 'gentrified_2016_2011TRUE')

res.lmm.0 <- update(res.lmm.zero, . ~ . - gentrified_2016_2011)
LRT(res.lmm.0, res.lmm.zero)

##     chi2_LR df p_value
## p_v       0  1       1

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.zero, plot = F)
  plot(simout, title = "DHARMa residual | Logit")
}

# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
  filter(area_esp_vert_2011ct != 0) %>%
  tidy_df(CT16, f)

res.lmm.cnt <- fitme(area_esp_vert_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: area_esp_vert_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                          Estimate Cond. SE t-value p-value
## (Intercept)               -0.8083  0.02176  -37.14       0
## gentrified_2016_2011TRUE  -0.4497  0.04445  -10.12       0
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##       1.rho 
## -0.01855287 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.2015  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -13407.45

#ci.lmm.cnt <- confint(res.lmm.cnt, 'gentrified_2016_2011TRUE')

res.lmm.0 <- update(res.lmm.cnt, . ~ . - gentrified_2016_2011)
LRT(res.lmm.0, res.lmm.cnt)

##      chi2_LR df p_value
## p_v -23.1309  1       1

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout, title = "DHARMa residual | Count (truncated Poisson)")
}

4.3.5 Canopy (trees)

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_high_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_high_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx, 
                 family = Poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_high_2011ct ~ gentrified_2016_2011 + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                          Estimate Cond. SE t-value  p-value
## (Intercept)               -1.6372  0.03951 -41.442 0.000000
## gentrified_2016_2011TRUE  -0.1505  0.04990  -3.015 0.002568
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1380114 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.1457  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2174.734

ci.lmm <- confint(res.lmm, 'gentrified_2016_2011TRUE')

## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE 
##                    -0.24862544                    -0.05240792

# Testing main effect
res.lmm.0 <- update(res.lmm, . ~ . - gentrified_2016_2011)
LRT(res.lmm.0, res.lmm)

##      chi2_LR df     p_value
## p_v 9.025009  1 0.002663107

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

Interpretation of the gentrified_2016_2011 coefficient (significant): we have \(e^{-0.15}\) = 0.86 (95%CI: 0.78 – 0.949), which means gentrified CTs exhibit a 14% decrease in the percentage of tree canopy area within CT in 2011

5 Association between SES and Urban Interventions at baseline (2011)

Looking at objective #1 | do urban interventions tend to be located in low SES neighborhoods?. We look at \[Urban Intervention_{2011 \to 2016} = f(SES_{2011})\] as well as \[Urban Intervention_{2011 \to 2016} = f(Gentrification_{2011 \to 2016})\]

Here \(Urban Intervention\) means the changes in the urban environment features, such as variation of bike lane length, greenness coverage, etc.

Apparently, change in rates are rather modeled as a pre-post model, as suggested in Jones’s answer. Hence, we model change in interventions between 2011 and 2016 as \[UC_{2016} = \beta_{0} + \beta_{1} * UC_{2011} + \beta_{2} * SES\]

Mostly redundant with multi-level approach below.

5.1 UI vs Pampalon | material

5.1.1 Bike lane length change

Bike lane length change

f <- Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                              Estimate Cond. SE t-value p-value
## (Intercept)                                 -2.927486 0.109112 -26.830 0.00000
## wSCOREMAT.2011.Q                            -0.068554 0.028786  -2.381 0.01724
## Bike_lane.by.street.2011ct                   0.065494 0.005882  11.135 0.00000
## wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct  0.003233 0.001796   1.801 0.07176
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1373553 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.29  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2271.208

ci.lmm.1 <- confint(res.lmm, "wSCOREMAT.2011.Q")

## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q 
##            -0.12531917            -0.01199128

ci.lmm.2 <- confint(res.lmm, "Bike_lane.by.street.2011ct")

## lower Bike_lane.by.street.2011ct upper Bike_lane.by.street.2011ct 
##                       0.05395184                       0.07726931

ci.lmm.3 <- confint(res.lmm, "wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct")

## lower wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct 
##                                     -0.0003019513 
## upper wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct 
##                                      0.0067687791

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_mat_bk <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_bk

# AIC
AIC(res.lmm)

Interpretation of the wSCOREMAT.2011.Q coefficient (significant): we have \(e^{-0.069}\) = 0.934 (95%CI: 0.882 – 0.988)

Coefficient interpretation (marginal + interaction) as relative risk. (NB need to be interpreted with a grain of salt given that interaction term is not significant)

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(wSCOREMAT.2011.Q=c(1,3,5))
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['Bike_lane.by.street.2011ct'] + i_e["wSCOREMAT.2011.Q"] * res.lmm[["fixef"]]['wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(Bike_lane.by.street.2011ct = quantile(clean_bei$df$Bike_lane.by.street.2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['wSCOREMAT.2011.Q'] + i_e["Bike_lane.by.street.2011ct"] * res.lmm[["fixef"]]['wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e

5.1.2 Bike lane length change (hurdle)

Bike lane length change

f <- Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm.0 <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = binomial())
summary(res.lmm.0, details = c(p_value="Wald"))

## formula: I(Bike_lane_total.hm.2016ct != 0) ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                                             Estimate Cond. SE  t-value
## (Intercept)                                 -4.84724   0.8762 -5.53209
## wSCOREMAT.2011.Q                            -0.02047   0.2051 -0.09978
## Bike_lane.by.street.2011ct                   3.65315   2.5340  1.44163
## wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct -0.31922   0.5908 -0.54034
##                                               p-value
## (Intercept)                                 3.164e-08
## wSCOREMAT.2011.Q                            9.205e-01
## Bike_lane.by.street.2011ct                  1.494e-01
## wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct 5.890e-01
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1377177 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  5.634  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -149.4448

# NB: not converging | but not needed, use Wald p-value instead or LRT
# ci.lmm.0.1 <- confint(res.lmm.0, "wSCOREMAT.2011.Q")
# ci.lmm.0.2 <- confint(res.lmm.0, "Bike_lane.by.street.2011ct")
# ci.lmm.0.3 <- confint(res.lmm.0, "wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
  plot(simout)
}

# Plotting interaction || MEANINGLESS
# .t <- clean_bei$df
# .t["predict"] <- as.vector(predict(res.lmm.0))
# 
# ggplot(.t, aes(y=predict, x=Bike_lane.by.street.2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) + 
#   geom_smooth(method = "lm", alpha=.1) +
#   geom_point() + 
#   labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')

# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
  filter(Bike_lane_total.hm.2016ct != 0) %>%
  tidy_df(CT16, f)

res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                               Estimate Cond. SE  t-value
## (Intercept)                                 -2.6497876 0.082688 -32.0455
## wSCOREMAT.2011.Q                            -0.0158610 0.022612  -0.7014
## Bike_lane.by.street.2011ct                   0.0527900 0.004345  12.1486
## wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct  0.0004125 0.001355   0.3045
##                                             p-value
## (Intercept)                                  0.0000
## wSCOREMAT.2011.Q                             0.4830
## Bike_lane.by.street.2011ct                   0.0000
## wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct  0.7608
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1502943 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.1123  
## # of obs: 564; # of groups: ct_no, 564 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1884.235

# ci.lmm.cnt.1 <- confint(res.lmm.cnt, "wSCOREMAT.2011.Q")
# ci.lmm.cnt.2 <- confint(res.lmm.cnt, "Bike_lane.by.street.2011ct")
# ci.lmm.cnt.3 <- confint(res.lmm.cnt, "wSCOREMAT.2011.Q:Bike_lane.by.street.2011ct")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei.cnt$df
.t["predict"] <- as.vector(predict(res.lmm.cnt))

g_mat_bk.hrdl <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1, linetype="dashed") +
  labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_bk.hrdl

# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)

# Combined AIC
"Combined Hurdle AIC:"

## [1] "Combined Hurdle AIC:"

-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)

## [1] 4079.36

5.1.3 Greenness change

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_2011ct + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_2017ct ~ wSCOREMAT.2011.Q*pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_2011ct + 
##     offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                       Estimate  Cond. SE t-value   p-value
## (Intercept)                          -1.458572 0.0484021 -30.134 0.0000000
## wSCOREMAT.2011.Q                     -0.057615 0.0141384  -4.075 0.0000460
## pct_esp_vert_2011ct                   0.016103 0.0008524  18.893 0.0000000
## wSCOREMAT.2011.Q:pct_esp_vert_2011ct  0.001066 0.0002922   3.648 0.0002641
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1363052 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.004542  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2010.443

ci.lmm.1 <- confint(res.lmm, "wSCOREMAT.2011.Q")

## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q 
##            -0.08539778            -0.02969523

ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_2011ct")

## lower pct_esp_vert_2011ct upper pct_esp_vert_2011ct 
##                0.01441588                0.01782357

ci.lmm.3 <- confint(res.lmm, "wSCOREMAT.2011.Q:pct_esp_vert_2011ct")

## lower wSCOREMAT.2011.Q:pct_esp_vert_2011ct 
##                               0.0004922499 
## upper wSCOREMAT.2011.Q:pct_esp_vert_2011ct 
##                               0.0016396399

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_mat_gr <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_gr

Interpretation of the wSCOREMAT.2011.Q coefficient (significant): we have \(e^{-0.058}\) = 0.944 (95%CI: 0.918 – 0.971)

Coefficient interpretation (marginal + interaction) as relative risk.

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(wSCOREMAT.2011.Q=c(1,3,5))
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_2011ct'] + i_e["wSCOREMAT.2011.Q"] * res.lmm[["fixef"]]['wSCOREMAT.2011.Q:pct_esp_vert_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_2011ct = quantile(clean_bei$df$pct_esp_vert_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['wSCOREMAT.2011.Q'] + i_e["pct_esp_vert_2011ct"] * res.lmm[["fixef"]]['wSCOREMAT.2011.Q:pct_esp_vert_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e

5.1.4 Canopy (trees) change

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_high_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_high_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_high_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_high_2011ct + 
##     offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                            Estimate  Cond. SE t-value   p-value
## (Intercept)                               -2.093050 0.0484735 -43.179 0.000e+00
## wSCOREMAT.2011.Q                          -0.114758 0.0151884  -7.556 4.174e-14
## pct_esp_vert_high_2011ct                   0.028402 0.0015954  17.802 0.000e+00
## wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct  0.005477 0.0006629   8.263 1.110e-16
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1379226 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.01035  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1790.562

ci.lmm.1 <- confint(res.lmm, "wSCOREMAT.2011.Q")

## lower wSCOREMAT.2011.Q upper wSCOREMAT.2011.Q 
##            -0.14459458            -0.08498978

ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_high_2011ct")

## lower pct_esp_vert_high_2011ct upper pct_esp_vert_high_2011ct 
##                     0.02527578                     0.03155832

ci.lmm.3 <- confint(res.lmm, "wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct")

## lower wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct 
##                                     0.004178959 
## upper wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct 
##                                     0.006787995

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_mat_tc <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_high_2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_tc

Interpretation of the wSCOREMAT.2011.Q coefficient (significant): we have \(e^{-0.115}\) = 0.892 (95%CI: 0.865 – 0.919)

Coefficient interpretation (marginal + interaction) as relative risk.

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(wSCOREMAT.2011.Q=c(1,3,5))
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_high_2011ct'] + i_e["wSCOREMAT.2011.Q"] * res.lmm[["fixef"]]['wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_high_2011ct = quantile(clean_bei$df$pct_esp_vert_high_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['wSCOREMAT.2011.Q'] + i_e["pct_esp_vert_high_2011ct"] * res.lmm[["fixef"]]['wSCOREMAT.2011.Q:pct_esp_vert_high_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e

5.2 UI vs Visible minority

5.2.1 Bike lane length change

Bike lane ratio to streets (in %). With or without controlling for UC in 2011 at the Census tract level

f <- Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                 Estimate Cond. SE t-value
## (Intercept)                                    -2.741242 0.147329 -18.606
## vis_minority_2011.Q                            -0.112340 0.038099  -2.949
## Bike_lane.by.street.2011ct                      0.053790 0.008465   6.354
## vis_minority_2011.Q:Bike_lane.by.street.2011ct  0.006389 0.002364   2.703
##                                                  p-value
## (Intercept)                                    0.000e+00
## vis_minority_2011.Q                            3.192e-03
## Bike_lane.by.street.2011ct                     2.094e-10
## vis_minority_2011.Q:Bike_lane.by.street.2011ct 6.882e-03
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1373916 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.2888  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2266.406

ci.lmm.1 <- confint(res.lmm, "vis_minority_2011.Q")

## lower vis_minority_2011.Q upper vis_minority_2011.Q 
##               -0.18768022               -0.03738233

ci.lmm.2 <- confint(res.lmm, "Bike_lane.by.street.2011ct")

## lower Bike_lane.by.street.2011ct upper Bike_lane.by.street.2011ct 
##                       0.03717189                       0.07052223

ci.lmm.3 <- confint(res.lmm, "vis_minority_2011.Q:Bike_lane.by.street.2011ct")

## lower vis_minority_2011.Q:Bike_lane.by.street.2011ct 
##                                          0.001745146 
## upper vis_minority_2011.Q:Bike_lane.by.street.2011ct 
##                                          0.011067714

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_vis_bk <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_bk

# AIC
AIC(res.lmm)

Interpretation of the vis_minority_2011.Q coefficient (significant): we have \(e^{-0.112}\) = 0.894 (95%CI: 0.829 – 0.963)

Coefficient interpretation (marginal + interaction) as relative risk.

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(vis_minority_2011.Q=c(1,3,5))
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['Bike_lane.by.street.2011ct'] + i_e["vis_minority_2011.Q"] * res.lmm[["fixef"]]['vis_minority_2011.Q:Bike_lane.by.street.2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(Bike_lane.by.street.2011ct = quantile(clean_bei$df$Bike_lane.by.street.2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['vis_minority_2011.Q'] + i_e["Bike_lane.by.street.2011ct"] * res.lmm[["fixef"]]['vis_minority_2011.Q:Bike_lane.by.street.2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e

5.2.2 Bike lane length change (hurdle)

Bike lane length change

f <- Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm.0 <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = binomial())
summary(res.lmm.0, details = c(p_value="Wald"))

## formula: I(Bike_lane_total.hm.2016ct != 0) ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                                                Estimate Cond. SE t-value
## (Intercept)                                     -2.8125   1.2554 -2.2403
## vis_minority_2011.Q                             -0.5283   0.2988 -1.7680
## Bike_lane.by.street.2011ct                       8.8040   7.8331  1.1240
## vis_minority_2011.Q:Bike_lane.by.street.2011ct  -1.3332   1.6546 -0.8057
##                                                p-value
## (Intercept)                                    0.02507
## vis_minority_2011.Q                            0.07705
## Bike_lane.by.street.2011ct                     0.26103
## vis_minority_2011.Q:Bike_lane.by.street.2011ct 0.42039
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1379107 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  6.666  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -145.2838

# NB: not converging
# ci.lmm.0.1 <- confint(res.lmm.0, "vis_minority_2011.Q")
# ci.lmm.0.2 <- confint(res.lmm.0, "Bike_lane.by.street.2011ct")
# ci.lmm.0.3 <- confint(res.lmm.0, "vis_minority_2011.Q:Bike_lane.by.street.2011ct")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
  plot(simout)
}

# Plotting interaction  || MEANINGLESS
# .t <- clean_bei$df
# .t["predict"] <- as.vector(predict(res.lmm.0))
# 
# ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) + 
#   geom_smooth(method = "lm", alpha=.1) +
#   labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')

# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
  filter(Bike_lane_total.hm.2016ct != 0) %>%
  tidy_df(CT16, f)

res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                  Estimate Cond. SE  t-value
## (Intercept)                                    -2.5785538 0.113994 -22.6201
## vis_minority_2011.Q                            -0.0332060 0.029423  -1.1286
## Bike_lane.by.street.2011ct                      0.0530570 0.006181   8.5842
## vis_minority_2011.Q:Bike_lane.by.street.2011ct  0.0002598 0.001755   0.1481
##                                                p-value
## (Intercept)                                     0.0000
## vis_minority_2011.Q                             0.2591
## Bike_lane.by.street.2011ct                      0.0000
## vis_minority_2011.Q:Bike_lane.by.street.2011ct  0.8823
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1504841 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.1118  
## # of obs: 564; # of groups: ct_no, 564 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1880.599

# ci.lmm.cnt.1 <- confint(res.lmm.cnt, "vis_minority_2011.Q")
# ci.lmm.cnt.2 <- confint(res.lmm.cnt, "Bike_lane.by.street.2011ct")
# ci.lmm.cnt.3 <- confint(res.lmm.cnt, "vis_minority_2011.Q:Bike_lane.by.street.2011ct")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei.cnt$df
.t["predict"] <- as.vector(predict(res.lmm.cnt))

g_vis_bk.hrdl <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1, linetype = "dashed") +
  labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_bk.hrdl

# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)

# Combined AIC
"Combined Hurdle AIC:"

## [1] "Combined Hurdle AIC:"

-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)

## [1] 4063.765

5.2.3 Greenness change

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_2017ct ~ vis_minority_2011.Q * pct_esp_vert_2011ct + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_2017ct ~ vis_minority_2011.Q * pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_2017ct ~ vis_minority_2011.Q * pct_esp_vert_2011ct + 
##     offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                          Estimate  Cond. SE t-value   p-value
## (Intercept)                             -1.330696 0.0651917 -20.412 0.000e+00
## vis_minority_2011.Q                     -0.085469 0.0174520  -4.897 9.711e-07
## pct_esp_vert_2011ct                      0.014220 0.0011258  12.631 0.000e+00
## vis_minority_2011.Q:pct_esp_vert_2011ct  0.001415 0.0003249   4.354 1.336e-05
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1360559 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.004362  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2006.168

ci.lmm.1 <- confint(res.lmm, "vis_minority_2011.Q")

## lower vis_minority_2011.Q upper vis_minority_2011.Q 
##               -0.11970466               -0.05096849

ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_2011ct")

## lower pct_esp_vert_2011ct upper pct_esp_vert_2011ct 
##                0.01200379                0.01646749

ci.lmm.3 <- confint(res.lmm, "vis_minority_2011.Q:pct_esp_vert_2011ct")

## lower vis_minority_2011.Q:pct_esp_vert_2011ct 
##                                  0.0007769562 
## upper vis_minority_2011.Q:pct_esp_vert_2011ct 
##                                  0.0020534515

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_vis_gr <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_gr

Interpretation of the vis_minority_2011.Q coefficient (significant): we have \(e^{-0.085}\) = 0.918 (95%CI: 0.887 – 0.95)

Coefficient interpretation (marginal + interaction) as relative risk.

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(vis_minority_2011.Q=c(1,3,5))
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_2011ct'] + i_e["vis_minority_2011.Q"] * res.lmm[["fixef"]]['vis_minority_2011.Q:pct_esp_vert_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_2011ct = quantile(clean_bei$df$pct_esp_vert_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['vis_minority_2011.Q'] + i_e["pct_esp_vert_2011ct"] * res.lmm[["fixef"]]['vis_minority_2011.Q:pct_esp_vert_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e

5.2.4 Canopy (trees) change

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_high_2017ct ~ vis_minority_2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_high_2017ct ~ vis_minority_2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_high_2017ct ~ vis_minority_2011.Q * pct_esp_vert_high_2011ct + 
##     offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                               Estimate  Cond. SE t-value
## (Intercept)                                  -2.038501 0.0682514 -29.868
## vis_minority_2011.Q                          -0.102183 0.0188781  -5.413
## pct_esp_vert_high_2011ct                      0.026803 0.0023166  11.570
## vis_minority_2011.Q:pct_esp_vert_high_2011ct  0.004335 0.0007174   6.043
##                                                p-value
## (Intercept)                                  0.000e+00
## vis_minority_2011.Q                          6.206e-08
## pct_esp_vert_high_2011ct                     0.000e+00
## vis_minority_2011.Q:pct_esp_vert_high_2011ct 1.517e-09
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1379804 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.01138  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1803.976

ci.lmm.1 <- confint(res.lmm, "vis_minority_2011.Q")

## lower vis_minority_2011.Q upper vis_minority_2011.Q 
##               -0.13931992               -0.06512398

ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_high_2011ct")

## lower pct_esp_vert_high_2011ct upper pct_esp_vert_high_2011ct 
##                     0.02225584                     0.03137141

ci.lmm.3 <- confint(res.lmm, "vis_minority_2011.Q:pct_esp_vert_high_2011ct")

## lower vis_minority_2011.Q:pct_esp_vert_high_2011ct 
##                                        0.002928509 
## upper vis_minority_2011.Q:pct_esp_vert_high_2011ct 
##                                        0.005753437

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_vis_tc <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_high_2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_tc

Interpretation of the vis_minority_2011.Q coefficient (significant): we have \(e^{-0.102}\) = 0.903 (95%CI: 0.87 – 0.937)

Coefficient interpretation (marginal + interaction) as relative risk.

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(vis_minority_2011.Q=c(1,3,5))
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_high_2011ct'] + i_e["vis_minority_2011.Q"] * res.lmm[["fixef"]]['vis_minority_2011.Q:pct_esp_vert_high_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_high_2011ct = quantile(clean_bei$df$pct_esp_vert_high_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['vis_minority_2011.Q'] + i_e["pct_esp_vert_high_2011ct"] * res.lmm[["fixef"]]['vis_minority_2011.Q:pct_esp_vert_high_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e

5.3 UI vs gentrified CT

Gentrified CT between 2011 and 2016

5.3.1 Bike lane length change

Bike lane ratio to streets (in %) at the Census tract level

f <- Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                      Estimate Cond. SE  t-value
## (Intercept)                                         -3.163664 0.067593 -46.8048
## gentrified_2016_2011TRUE                             0.076678 0.089229   0.8593
## Bike_lane.by.street.2011ct                           0.077872 0.003409  22.8434
## gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct -0.006855 0.005835  -1.1749
##                                                     p-value
## (Intercept)                                          0.0000
## gentrified_2016_2011TRUE                             0.3902
## Bike_lane.by.street.2011ct                           0.0000
## gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct  0.2400
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1375707 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.2922  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2270.491

ci.lmm.1 <- confint(res.lmm, "gentrified_2016_2011TRUE")

## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE 
##                     -0.1000857                      0.2518371

ci.lmm.2 <- confint(res.lmm, "Bike_lane.by.street.2011ct")

## lower Bike_lane.by.street.2011ct upper Bike_lane.by.street.2011ct 
##                       0.07111801                       0.08484937

ci.lmm.3 <- confint(res.lmm, "gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct")

## lower gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct 
##                                              -0.018320553 
## upper gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct 
##                                               0.004656173

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_gen_bk <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(gentrified_2016_2011), fill=factor(gentrified_2016_2011))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='gentrified_2016_2011', fill='gentrified_2016_2011')
g_gen_bk

# AIC
AIC(res.lmm)

Interpretation of the gentrified_2016_2011 coefficient (significant): we have \(e^{0.077}\) = 1.08 (95%CI: 0.905 – 1.286)

Coefficient interpretation (marginal + interaction) as relative risk. (NB need to be interpreted with a grain of salt given that gentrified term and interaction term are not significant)

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(gentrified_2016_2011=0:1)
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['Bike_lane.by.street.2011ct'] + i_e["gentrified_2016_2011"] * res.lmm[["fixef"]]['gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(Bike_lane.by.street.2011ct = quantile(clean_bei$df$Bike_lane.by.street.2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['gentrified_2016_2011TRUE'] + i_e["Bike_lane.by.street.2011ct"] * res.lmm[["fixef"]]['gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e

5.3.2 Bike lane length change (hurdle)

Bike lane length change

f <- Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm.0 <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = binomial())
summary(res.lmm.0, details = c(p_value="Wald"))

## formula: I(Bike_lane_total.hm.2016ct != 0) ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                                                     Estimate Cond. SE t-value
## (Intercept)                                          -4.9123   0.5013 -9.7998
## gentrified_2016_2011TRUE                              0.1917   0.5944  0.3226
## Bike_lane.by.street.2011ct                            3.4688   1.4200  2.4427
## gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct  -1.6031   1.6952 -0.9457
##                                                     p-value
## (Intercept)                                         0.00000
## gentrified_2016_2011TRUE                            0.74703
## Bike_lane.by.street.2011ct                          0.01458
## gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct 0.34432
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1379976 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  6.356  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -146.8066

# NB: not converging
# ci.lmm.0.1 <- confint(res.lmm.0, "gentrified_2016_2011TRUE")
# ci.lmm.0.2 <- confint(res.lmm.0, "Bike_lane.by.street.2011ct")
# ci.lmm.0.3 <- confint(res.lmm.0, "gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
  plot(simout)
}

# Plotting interaction || MEANINGLESS
# .t <- clean_bei$df
# .t["predict"] <- as.vector(predict(res.lmm.0))
# 
# ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(gentrified_2016_2011), fill=factor(gentrified_2016_2011))) + 
#   geom_smooth(method = "lm", alpha=.1) +
#   labs(colour='gentrified_2016_2011', fill='gentrified_2016_2011')

# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
  filter(Bike_lane_total.hm.2016ct != 0) %>%
  tidy_df(CT16, f)

res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2016ct ~ gentrified_2016_2011 * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                     Estimate Cond. SE t-value
## (Intercept)                                         -2.75518 0.046857 -58.800
## gentrified_2016_2011TRUE                             0.20178 0.071981   2.803
## Bike_lane.by.street.2011ct                           0.05753 0.002441  23.564
## gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct -0.01157 0.004373  -2.646
##                                                      p-value
## (Intercept)                                         0.000000
## gentrified_2016_2011TRUE                            0.005060
## Bike_lane.by.street.2011ct                          0.000000
## gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct 0.008136
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1504352 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.1088  
## # of obs: 564; # of groups: ct_no, 564 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1877.541

# ci.lmm.cnt.1 <- confint(res.lmm.cnt, "gentrified_2016_2011TRUE")
# ci.lmm.cnt.2 <- confint(res.lmm.cnt, "Bike_lane.by.street.2011ct")
# ci.lmm.cnt.3 <- confint(res.lmm.cnt, "gentrified_2016_2011TRUE:Bike_lane.by.street.2011ct")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei.cnt$df
.t["predict"] <- as.vector(predict(res.lmm.cnt))

g_gen_bk.hrdl <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(gentrified_2016_2011), fill=factor(gentrified_2016_2011))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='gentrified_2016_2011', fill='gentrified_2016_2011')
g_gen_bk.hrdl

# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)

# Combined AIC
"Combined Hurdle AIC:"

## [1] "Combined Hurdle AIC:"

-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)

## [1] 4060.696

5.3.3 Greenness change

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_2017ct ~ gentrified_2016_2011 * pct_esp_vert_2011ct + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_2017ct ~ gentrified_2016_2011 * pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_2017ct ~ gentrified_2016_2011 * pct_esp_vert_2011ct + 
##     offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                               Estimate  Cond. SE t-value
## (Intercept)                                  -1.613014 0.0257795 -62.570
## gentrified_2016_2011TRUE                     -0.202554 0.0643593  -3.147
## pct_esp_vert_2011ct                           0.018742 0.0004559  41.105
## gentrified_2016_2011TRUE:pct_esp_vert_2011ct  0.005263 0.0016871   3.120
##                                               p-value
## (Intercept)                                  0.000000
## gentrified_2016_2011TRUE                     0.001648
## pct_esp_vert_2011ct                          0.000000
## gentrified_2016_2011TRUE:pct_esp_vert_2011ct 0.001811
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1349448 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.004838  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2012.743

ci.lmm.1 <- confint(res.lmm, "gentrified_2016_2011TRUE")

## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE 
##                    -0.32924546                    -0.07617365

ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_2011ct")

## lower pct_esp_vert_2011ct upper pct_esp_vert_2011ct 
##                0.01782698                0.01968221

ci.lmm.3 <- confint(res.lmm, "gentrified_2016_2011TRUE:pct_esp_vert_2011ct")

## lower gentrified_2016_2011TRUE:pct_esp_vert_2011ct 
##                                        0.001943885 
## upper gentrified_2016_2011TRUE:pct_esp_vert_2011ct 
##                                        0.008578808

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_gen_gr <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_2011ct, colour=factor(gentrified_2016_2011), fill=factor(gentrified_2016_2011))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='gentrified_2016_2011', fill='gentrified_2016_2011')
g_gen_gr

Interpretation of the gentrified_2016_2011 coefficient (significant): we have \(e^{-0.203}\) = 0.817 (95%CI: 0.719 – 0.927)

Coefficient interpretation (marginal + interaction) as relative risk.

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(gentrified_2016_2011=0:1)
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_2011ct'] + i_e["gentrified_2016_2011"] * res.lmm[["fixef"]]['gentrified_2016_2011TRUE:pct_esp_vert_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_2011ct = quantile(clean_bei$df$pct_esp_vert_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['gentrified_2016_2011TRUE'] + i_e["pct_esp_vert_2011ct"] * res.lmm[["fixef"]]['gentrified_2016_2011TRUE:pct_esp_vert_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e

5.3.4 Canopy (trees) change

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_high_2017ct ~ gentrified_2016_2011 * pct_esp_vert_high_2011ct + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_high_2017ct ~ gentrified_2016_2011 * pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_high_2017ct ~ gentrified_2016_2011 * pct_esp_vert_high_2011ct + 
##     offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                   Estimate  Cond. SE t-value
## (Intercept)                                       -2.35730 0.0279340 -84.388
## gentrified_2016_2011TRUE                          -0.20459 0.0715883  -2.858
## pct_esp_vert_high_2011ct                           0.03939 0.0009182  42.893
## gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct  0.01322 0.0035508   3.722
##                                                     p-value
## (Intercept)                                       0.0000000
## gentrified_2016_2011TRUE                          0.0042654
## pct_esp_vert_high_2011ct                          0.0000000
## gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct 0.0001974
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1380802 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.01222  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1814.314

ci.lmm.1 <- confint(res.lmm, "gentrified_2016_2011TRUE")

## lower gentrified_2016_2011TRUE upper gentrified_2016_2011TRUE 
##                    -0.34539306                    -0.06456076

ci.lmm.2 <- confint(res.lmm, "pct_esp_vert_high_2011ct")

## lower pct_esp_vert_high_2011ct upper pct_esp_vert_high_2011ct 
##                     0.03756862                     0.04126454

ci.lmm.3 <- confint(res.lmm, "gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct")

## lower gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct 
##                                             0.006232907 
## upper gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct 
##                                             0.020169474

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_gen_tc <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_high_2011ct, colour=factor(gentrified_2016_2011), fill=factor(gentrified_2016_2011))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='gentrified_2016_2011', fill='gentrified_2016_2011')
g_gen_tc

Interpretation of the gentrified_2016_2011 coefficient (significant): we have \(e^{-0.205}\) = 0.815 (95%CI: 0.708 – 0.937)

Coefficient interpretation (marginal + interaction) as relative risk.

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(gentrified_2016_2011=0:1)
i_e["UC2011.effect"] <- res.lmm[["fixef"]]['pct_esp_vert_high_2011ct'] + i_e["gentrified_2016_2011"] * res.lmm[["fixef"]]['gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct']
i_e["UC2011.relative.risk"] <- exp(i_e["UC2011.effect"])
i_e

# Compute combined effects (marginal + interaction) for each quintile of SES:
i_e <- data.frame(pct_esp_vert_high_2011ct = quantile(clean_bei$df$pct_esp_vert_high_2011ct)[2:4])
i_e["SES.effect"] <- res.lmm[["fixef"]]['gentrified_2016_2011TRUE'] + i_e["pct_esp_vert_high_2011ct"] * res.lmm[["fixef"]]['gentrified_2016_2011TRUE:pct_esp_vert_high_2011ct']
i_e["SES.relative.risk"] <- exp(i_e["SES.effect"])
i_e

6 Association between SES and UI | SES as categories

Following Yan and Behzad’s request, the exact same models as above are run, this time using SES quintiles as categories instead of ordinal values.

6.1 UI vs Pampalon | material

6.1.1 Bike lane length change

Bike lane length change

f <- Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(Bike_lane_total.hm.2016ct ~ factor(wSCOREMAT.2011.Q) * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2016ct ~ factor(wSCOREMAT.2011.Q) * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                        Estimate Cond. SE
## (Intercept)                                          -3.0197649 0.101886
## factor(wSCOREMAT.2011.Q)2                            -0.0156279 0.124989
## factor(wSCOREMAT.2011.Q)3                            -0.1832813 0.134330
## factor(wSCOREMAT.2011.Q)4                            -0.0595746 0.129557
## factor(wSCOREMAT.2011.Q)5                            -0.2994914 0.124780
## Bike_lane.by.street.2011ct                            0.0697205 0.005064
## factor(wSCOREMAT.2011.Q)2:Bike_lane.by.street.2011ct -0.0018757 0.007879
## factor(wSCOREMAT.2011.Q)3:Bike_lane.by.street.2011ct  0.0179516 0.008873
## factor(wSCOREMAT.2011.Q)4:Bike_lane.by.street.2011ct -0.0008867 0.007882
## factor(wSCOREMAT.2011.Q)5:Bike_lane.by.street.2011ct  0.0148520 0.007874
##                                                       t-value p-value
## (Intercept)                                          -29.6388 0.00000
## factor(wSCOREMAT.2011.Q)2                             -0.1250 0.90050
## factor(wSCOREMAT.2011.Q)3                             -1.3644 0.17244
## factor(wSCOREMAT.2011.Q)4                             -0.4598 0.64564
## factor(wSCOREMAT.2011.Q)5                             -2.4001 0.01639
## Bike_lane.by.street.2011ct                            13.7680 0.00000
## factor(wSCOREMAT.2011.Q)2:Bike_lane.by.street.2011ct  -0.2381 0.81183
## factor(wSCOREMAT.2011.Q)3:Bike_lane.by.street.2011ct   2.0232 0.04305
## factor(wSCOREMAT.2011.Q)4:Bike_lane.by.street.2011ct  -0.1125 0.91043
## factor(wSCOREMAT.2011.Q)5:Bike_lane.by.street.2011ct   1.8861 0.05928
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1373537 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.2855  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2267.644

# for (i in 2:5) {
#   ci.lmm.1 <- confint(res.lmm, paste0("factor(wSCOREMAT.2011.Q)", i))
# }
# ci.lmm.2 <- confint(res.lmm, "Bike_lane.by.street.2011ct")
# for (i in 2:5) {
#   ci.lmm.3 <- confint(res.lmm, paste0("factor(wSCOREMAT.2011.Q)", i, ":Bike_lane.by.street.2011ct"))
# }

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_mat_bk <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_bk

# AIC
AIC(res.lmm)

6.1.2 Bike lane length change (hurdle)

Bike lane length change

f <- Bike_lane_total.hm.2016ct ~ wSCOREMAT.2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm.0 <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ factor(wSCOREMAT.2011.Q) * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = binomial())
summary(res.lmm.0, details = c(p_value="Wald"))

## formula: I(Bike_lane_total.hm.2016ct != 0) ~ factor(wSCOREMAT.2011.Q) * 
##     Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                                                      Estimate Cond. SE t-value
## (Intercept)                                           -6.1432   0.8734 -7.0336
## factor(wSCOREMAT.2011.Q)2                              1.8520   0.9585  1.9322
## factor(wSCOREMAT.2011.Q)3                              2.2186   1.0634  2.0864
## factor(wSCOREMAT.2011.Q)4                              1.7289   0.9782  1.7673
## factor(wSCOREMAT.2011.Q)5                              0.7182   0.9227  0.7784
## Bike_lane.by.street.2011ct                            51.7964 127.9150  0.4049
## factor(wSCOREMAT.2011.Q)2:Bike_lane.by.street.2011ct -50.3057 127.9176 -0.3933
## factor(wSCOREMAT.2011.Q)3:Bike_lane.by.street.2011ct -47.5550 128.2412 -0.3708
## factor(wSCOREMAT.2011.Q)4:Bike_lane.by.street.2011ct 159.9890 891.0603  0.1795
## factor(wSCOREMAT.2011.Q)5:Bike_lane.by.street.2011ct -50.0434 127.9165 -0.3912
##                                                        p-value
## (Intercept)                                          2.013e-12
## factor(wSCOREMAT.2011.Q)2                            5.333e-02
## factor(wSCOREMAT.2011.Q)3                            3.694e-02
## factor(wSCOREMAT.2011.Q)4                            7.718e-02
## factor(wSCOREMAT.2011.Q)5                            4.363e-01
## Bike_lane.by.street.2011ct                           6.855e-01
## factor(wSCOREMAT.2011.Q)2:Bike_lane.by.street.2011ct 6.941e-01
## factor(wSCOREMAT.2011.Q)3:Bike_lane.by.street.2011ct 7.108e-01
## factor(wSCOREMAT.2011.Q)4:Bike_lane.by.street.2011ct 8.575e-01
## factor(wSCOREMAT.2011.Q)5:Bike_lane.by.street.2011ct 6.956e-01
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1376862 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  4.612  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -142.7166

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
  plot(simout)
}

# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
  filter(Bike_lane_total.hm.2016ct != 0) %>%
  tidy_df(CT16, f)

res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ factor(wSCOREMAT.2011.Q) * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2016ct ~ factor(wSCOREMAT.2011.Q) * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                        Estimate Cond. SE
## (Intercept)                                          -2.6145832 0.074587
## factor(wSCOREMAT.2011.Q)2                            -0.0480572 0.096204
## factor(wSCOREMAT.2011.Q)3                            -0.2912452 0.101969
## factor(wSCOREMAT.2011.Q)4                            -0.1329661 0.098916
## factor(wSCOREMAT.2011.Q)5                            -0.0300572 0.097780
## Bike_lane.by.street.2011ct                            0.0506292 0.003668
## factor(wSCOREMAT.2011.Q)2:Bike_lane.by.street.2011ct  0.0001294 0.005675
## factor(wSCOREMAT.2011.Q)3:Bike_lane.by.street.2011ct  0.0212050 0.006469
## factor(wSCOREMAT.2011.Q)4:Bike_lane.by.street.2011ct  0.0044888 0.005729
## factor(wSCOREMAT.2011.Q)5:Bike_lane.by.street.2011ct -0.0019039 0.006090
##                                                       t-value  p-value
## (Intercept)                                          -35.0543 0.000000
## factor(wSCOREMAT.2011.Q)2                             -0.4995 0.617404
## factor(wSCOREMAT.2011.Q)3                             -2.8562 0.004287
## factor(wSCOREMAT.2011.Q)4                             -1.3442 0.178875
## factor(wSCOREMAT.2011.Q)5                             -0.3074 0.758543
## Bike_lane.by.street.2011ct                            13.8020 0.000000
## factor(wSCOREMAT.2011.Q)2:Bike_lane.by.street.2011ct   0.0228 0.981812
## factor(wSCOREMAT.2011.Q)3:Bike_lane.by.street.2011ct   3.2781 0.001045
## factor(wSCOREMAT.2011.Q)4:Bike_lane.by.street.2011ct   0.7835 0.433317
## factor(wSCOREMAT.2011.Q)5:Bike_lane.by.street.2011ct  -0.3126 0.754587
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1503544 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.1079  
## # of obs: 564; # of groups: ct_no, 564 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1877.346

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei.cnt$df
.t["predict"] <- as.vector(predict(res.lmm.cnt))

g_mat_bk.hrdl <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1, linetype="dashed") +
  labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_bk.hrdl

# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)

# Combined AIC
"Combined Hurdle AIC:"

## [1] "Combined Hurdle AIC:"

-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)

## [1] 4064.124

6.1.3 Greenness change

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_2011ct + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_2017ct ~ factor(wSCOREMAT.2011.Q)* pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_2017ct ~ factor(wSCOREMAT.2011.Q) * pct_esp_vert_2011ct + 
##     offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                 Estimate  Cond. SE  t-value
## (Intercept)                                   -1.5072099 0.0449681 -33.5173
## factor(wSCOREMAT.2011.Q)2                     -0.1205581 0.0634337  -1.9005
## factor(wSCOREMAT.2011.Q)3                     -0.0708225 0.0649834  -1.0899
## factor(wSCOREMAT.2011.Q)4                     -0.1380149 0.0719175  -1.9191
## factor(wSCOREMAT.2011.Q)5                     -0.3096978 0.0657409  -4.7109
## pct_esp_vert_2011ct                            0.0172320 0.0007454  23.1175
## factor(wSCOREMAT.2011.Q)2:pct_esp_vert_2011ct  0.0015789 0.0010825   1.4586
## factor(wSCOREMAT.2011.Q)3:pct_esp_vert_2011ct  0.0008741 0.0011405   0.7664
## factor(wSCOREMAT.2011.Q)4:pct_esp_vert_2011ct  0.0024727 0.0014599   1.6938
## factor(wSCOREMAT.2011.Q)5:pct_esp_vert_2011ct  0.0061844 0.0014984   4.1274
##                                                 p-value
## (Intercept)                                   0.000e+00
## factor(wSCOREMAT.2011.Q)2                     5.736e-02
## factor(wSCOREMAT.2011.Q)3                     2.758e-01
## factor(wSCOREMAT.2011.Q)4                     5.498e-02
## factor(wSCOREMAT.2011.Q)5                     2.466e-06
## pct_esp_vert_2011ct                           0.000e+00
## factor(wSCOREMAT.2011.Q)2:pct_esp_vert_2011ct 1.447e-01
## factor(wSCOREMAT.2011.Q)3:pct_esp_vert_2011ct 4.434e-01
## factor(wSCOREMAT.2011.Q)4:pct_esp_vert_2011ct 9.030e-02
## factor(wSCOREMAT.2011.Q)5:pct_esp_vert_2011ct 3.669e-05
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1365155 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.004416  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2005.519

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_mat_gr <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_gr

6.1.4 Canopy (trees) change

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_high_2017ct ~ wSCOREMAT.2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_high_2017ct ~ factor(wSCOREMAT.2011.Q) * pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_high_2017ct ~ factor(wSCOREMAT.2011.Q) * pct_esp_vert_high_2011ct + 
##     offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                     Estimate Cond. SE t-value
## (Intercept)                                        -2.205007 0.041604 -53.000
## factor(wSCOREMAT.2011.Q)2                          -0.170448 0.063635  -2.679
## factor(wSCOREMAT.2011.Q)3                          -0.130104 0.066325  -1.962
## factor(wSCOREMAT.2011.Q)4                          -0.283483 0.071434  -3.968
## factor(wSCOREMAT.2011.Q)5                          -0.694646 0.077043  -9.016
## pct_esp_vert_high_2011ct                            0.034056 0.001196  28.474
## factor(wSCOREMAT.2011.Q)2:pct_esp_vert_high_2011ct  0.007114 0.002141   3.322
## factor(wSCOREMAT.2011.Q)3:pct_esp_vert_high_2011ct  0.004798 0.002391   2.007
## factor(wSCOREMAT.2011.Q)4:pct_esp_vert_high_2011ct  0.013497 0.003052   4.422
## factor(wSCOREMAT.2011.Q)5:pct_esp_vert_high_2011ct  0.038485 0.004293   8.964
##                                                      p-value
## (Intercept)                                        0.000e+00
## factor(wSCOREMAT.2011.Q)2                          7.395e-03
## factor(wSCOREMAT.2011.Q)3                          4.981e-02
## factor(wSCOREMAT.2011.Q)4                          7.233e-05
## factor(wSCOREMAT.2011.Q)5                          0.000e+00
## pct_esp_vert_high_2011ct                           0.000e+00
## factor(wSCOREMAT.2011.Q)2:pct_esp_vert_high_2011ct 8.941e-04
## factor(wSCOREMAT.2011.Q)3:pct_esp_vert_high_2011ct 4.475e-02
## factor(wSCOREMAT.2011.Q)4:pct_esp_vert_high_2011ct 9.769e-06
## factor(wSCOREMAT.2011.Q)5:pct_esp_vert_high_2011ct 0.000e+00
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1380144 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.00891  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1776.146

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_mat_tc <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_high_2011ct, colour=factor(wSCOREMAT.2011.Q), fill=factor(wSCOREMAT.2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='wSCOREMAT.2011.Q', fill='wSCOREMAT.2011.Q')
g_mat_tc

6.2 UI vs Visible minority

6.2.1 Bike lane length change

Bike lane ratio to streets (in %). With or without controlling for UC in 2011 at the Census tract level

f <- Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(Bike_lane_total.hm.2016ct ~ factor(vis_minority_2011.Q) * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2016ct ~ factor(vis_minority_2011.Q) * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                         Estimate Cond. SE
## (Intercept)                                             -2.41522  0.23066
## factor(vis_minority_2011.Q)2                            -0.66165  0.24390
## factor(vis_minority_2011.Q)3                            -0.62881  0.24137
## factor(vis_minority_2011.Q)4                            -0.86128  0.24335
## factor(vis_minority_2011.Q)5                            -0.84751  0.24407
## Bike_lane.by.street.2011ct                               0.04661  0.01223
## factor(vis_minority_2011.Q)2:Bike_lane.by.street.2011ct  0.02281  0.01347
## factor(vis_minority_2011.Q)3:Bike_lane.by.street.2011ct  0.02192  0.01317
## factor(vis_minority_2011.Q)4:Bike_lane.by.street.2011ct  0.04345  0.01363
## factor(vis_minority_2011.Q)5:Bike_lane.by.street.2011ct  0.03408  0.01351
##                                                         t-value   p-value
## (Intercept)                                             -10.471 0.0000000
## factor(vis_minority_2011.Q)2                             -2.713 0.0066723
## factor(vis_minority_2011.Q)3                             -2.605 0.0091822
## factor(vis_minority_2011.Q)4                             -3.539 0.0004013
## factor(vis_minority_2011.Q)5                             -3.472 0.0005158
## Bike_lane.by.street.2011ct                                3.811 0.0001385
## factor(vis_minority_2011.Q)2:Bike_lane.by.street.2011ct   1.694 0.0902987
## factor(vis_minority_2011.Q)3:Bike_lane.by.street.2011ct   1.664 0.0961098
## factor(vis_minority_2011.Q)4:Bike_lane.by.street.2011ct   3.189 0.0014275
## factor(vis_minority_2011.Q)5:Bike_lane.by.street.2011ct   2.522 0.0116824
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1374263 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.2808  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -2260.825

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_vis_bk <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_bk

# AIC
AIC(res.lmm)

6.2.2 Bike lane length change (hurdle)

Bike lane length change

f <- Bike_lane_total.hm.2016ct ~ vis_minority_2011.Q * Bike_lane.by.street.2011ct + offset(log(street_length.hm))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm.0 <- fitme(I(Bike_lane_total.hm.2016ct != 0) ~ factor(vis_minority_2011.Q) * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = binomial())
summary(res.lmm.0, details = c(p_value="Wald"))

## formula: I(Bike_lane_total.hm.2016ct != 0) ~ factor(vis_minority_2011.Q) * 
##     Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                                                         Estimate Cond. SE
## (Intercept)                                              -3.7597    2.418
## factor(vis_minority_2011.Q)2                             -1.3036    2.349
## factor(vis_minority_2011.Q)3                             -0.1278    2.442
## factor(vis_minority_2011.Q)4                             -0.9180    2.431
## factor(vis_minority_2011.Q)5                             -1.9497    2.440
## Bike_lane.by.street.2011ct                               91.8848 7989.814
## factor(vis_minority_2011.Q)2:Bike_lane.by.street.2011ct -86.9785 7989.916
## factor(vis_minority_2011.Q)3:Bike_lane.by.street.2011ct 133.7796 8135.766
## factor(vis_minority_2011.Q)4:Bike_lane.by.street.2011ct -88.9671 7989.814
## factor(vis_minority_2011.Q)5:Bike_lane.by.street.2011ct -89.7762 7989.814
##                                                          t-value p-value
## (Intercept)                                             -1.55491  0.1200
## factor(vis_minority_2011.Q)2                            -0.55500  0.5789
## factor(vis_minority_2011.Q)3                            -0.05232  0.9583
## factor(vis_minority_2011.Q)4                            -0.37761  0.7057
## factor(vis_minority_2011.Q)5                            -0.79914  0.4242
## Bike_lane.by.street.2011ct                               0.01150  0.9908
## factor(vis_minority_2011.Q)2:Bike_lane.by.street.2011ct -0.01089  0.9913
## factor(vis_minority_2011.Q)3:Bike_lane.by.street.2011ct  0.01644  0.9869
## factor(vis_minority_2011.Q)4:Bike_lane.by.street.2011ct -0.01114  0.9911
## factor(vis_minority_2011.Q)5:Bike_lane.by.street.2011ct -0.01124  0.9910
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1379579 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  5.382  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                        logLik
## p_v(h) (marginal L): -143.109

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
  plot(simout)
}

# Count model (Poisson)
clean_bei.cnt <- bei_df_aoi %>%
  filter(Bike_lane_total.hm.2016ct != 0) %>%
  tidy_df(CT16, f)

res.lmm.cnt <- fitme(Bike_lane_total.hm.2016ct ~ factor(vis_minority_2011.Q) * Bike_lane.by.street.2011ct + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = clean_bei.cnt$df, adjMatrix = clean_bei.cnt$nbmatx,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm.2016ct ~ factor(vis_minority_2011.Q) * Bike_lane.by.street.2011ct + 
##     offset(log(street_length.hm)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                          Estimate Cond. SE
## (Intercept)                                             -2.272331 0.164239
## factor(vis_minority_2011.Q)2                            -0.391292 0.178544
## factor(vis_minority_2011.Q)3                            -0.461271 0.173976
## factor(vis_minority_2011.Q)4                            -0.561317 0.176619
## factor(vis_minority_2011.Q)5                            -0.384225 0.176549
## Bike_lane.by.street.2011ct                               0.042063 0.008256
## factor(vis_minority_2011.Q)2:Bike_lane.by.street.2011ct  0.009861 0.009288
## factor(vis_minority_2011.Q)3:Bike_lane.by.street.2011ct  0.013153 0.008959
## factor(vis_minority_2011.Q)4:Bike_lane.by.street.2011ct  0.022163 0.009423
## factor(vis_minority_2011.Q)5:Bike_lane.by.street.2011ct  0.005333 0.009436
##                                                          t-value   p-value
## (Intercept)                                             -13.8355 0.000e+00
## factor(vis_minority_2011.Q)2                             -2.1916 2.841e-02
## factor(vis_minority_2011.Q)3                             -2.6513 8.017e-03
## factor(vis_minority_2011.Q)4                             -3.1781 1.482e-03
## factor(vis_minority_2011.Q)5                             -2.1763 2.953e-02
## Bike_lane.by.street.2011ct                                5.0950 3.488e-07
## factor(vis_minority_2011.Q)2:Bike_lane.by.street.2011ct   1.0617 2.884e-01
## factor(vis_minority_2011.Q)3:Bike_lane.by.street.2011ct   1.4680 1.421e-01
## factor(vis_minority_2011.Q)4:Bike_lane.by.street.2011ct   2.3521 1.867e-02
## factor(vis_minority_2011.Q)5:Bike_lane.by.street.2011ct   0.5652 5.719e-01
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1505247 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.1069  
## # of obs: 564; # of groups: ct_no, 564 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1873.797

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei.cnt$df
.t["predict"] <- as.vector(predict(res.lmm.cnt))

g_vis_bk.hrdl <- ggplot(.t, aes(y=100*predict/street_length.hm, x=Bike_lane.by.street.2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1, linetype = "dashed") +
  labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_bk.hrdl

# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)

# Combined AIC
"Combined Hurdle AIC:"

## [1] "Combined Hurdle AIC:"

-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)

## [1] 4057.812

6.2.3 Greenness change

Measuring canopy (i.e. greenness ~ grass & trees) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_2017ct ~ vis_minority_2011.Q * pct_esp_vert_2011ct + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_2017ct ~ factor(vis_minority_2011.Q) * pct_esp_vert_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_2017ct ~ factor(vis_minority_2011.Q) * pct_esp_vert_2011ct + 
##     offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                   Estimate Cond. SE t-value
## (Intercept)                                      -0.957726 0.095779  -9.999
## factor(vis_minority_2011.Q)2                     -0.637454 0.107212  -5.946
## factor(vis_minority_2011.Q)3                     -0.642669 0.102819  -6.251
## factor(vis_minority_2011.Q)4                     -0.697627 0.104089  -6.702
## factor(vis_minority_2011.Q)5                     -0.811202 0.104033  -7.798
## pct_esp_vert_2011ct                               0.009706 0.001507   6.439
## factor(vis_minority_2011.Q)2:pct_esp_vert_2011ct  0.008455 0.001714   4.932
## factor(vis_minority_2011.Q)3:pct_esp_vert_2011ct  0.008674 0.001641   5.285
## factor(vis_minority_2011.Q)4:pct_esp_vert_2011ct  0.009757 0.001721   5.669
## factor(vis_minority_2011.Q)5:pct_esp_vert_2011ct  0.012089 0.001762   6.860
##                                                    p-value
## (Intercept)                                      0.000e+00
## factor(vis_minority_2011.Q)2                     2.753e-09
## factor(vis_minority_2011.Q)3                     4.091e-10
## factor(vis_minority_2011.Q)4                     2.053e-11
## factor(vis_minority_2011.Q)5                     6.328e-15
## pct_esp_vert_2011ct                              1.204e-10
## factor(vis_minority_2011.Q)2:pct_esp_vert_2011ct 8.122e-07
## factor(vis_minority_2011.Q)3:pct_esp_vert_2011ct 1.255e-07
## factor(vis_minority_2011.Q)4:pct_esp_vert_2011ct 1.440e-08
## factor(vis_minority_2011.Q)5:pct_esp_vert_2011ct 6.909e-12
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1311812 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.003348  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1988.675

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_vis_gr <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_gr

6.2.4 Canopy (trees) change

Measuring high canopy (i.e. trees only) ratio within CT/buffer in 2011 (in %) at the Census tract level

f <- area_esp_vert_high_2017ct ~ vis_minority_2011.Q * pct_esp_vert_high_2011ct + offset(log(ct_area_hct))

clean_bei <- tidy_df(bei_df_aoi, CT16, f)

res.lmm <- fitme(area_esp_vert_high_2017ct ~ factor(vis_minority_2011.Q) * pct_esp_vert_high_2011ct + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = clean_bei$df, adjMatrix = clean_bei$nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_high_2017ct ~ factor(vis_minority_2011.Q) * pct_esp_vert_high_2011ct + 
##     offset(log(ct_area_hct)) + adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                                                         Estimate Cond. SE
## (Intercept)                                           -2.240e+00 0.081603
## factor(vis_minority_2011.Q)2                          -1.306e-02 0.095564
## factor(vis_minority_2011.Q)3                          -3.859e-02 0.091173
## factor(vis_minority_2011.Q)4                          -2.393e-01 0.094751
## factor(vis_minority_2011.Q)5                          -3.728e-01 0.095593
## pct_esp_vert_high_2011ct                               3.568e-02 0.002635
## factor(vis_minority_2011.Q)2:pct_esp_vert_high_2011ct  9.271e-05 0.003113
## factor(vis_minority_2011.Q)3:pct_esp_vert_high_2011ct  6.041e-04 0.002996
## factor(vis_minority_2011.Q)4:pct_esp_vert_high_2011ct  1.047e-02 0.003385
## factor(vis_minority_2011.Q)5:pct_esp_vert_high_2011ct  1.608e-02 0.003511
##                                                         t-value   p-value
## (Intercept)                                           -27.45366 0.000e+00
## factor(vis_minority_2011.Q)2                           -0.13661 8.913e-01
## factor(vis_minority_2011.Q)3                           -0.42331 6.721e-01
## factor(vis_minority_2011.Q)4                           -2.52536 1.156e-02
## factor(vis_minority_2011.Q)5                           -3.90005 9.617e-05
## pct_esp_vert_high_2011ct                               13.54153 0.000e+00
## factor(vis_minority_2011.Q)2:pct_esp_vert_high_2011ct   0.02978 9.762e-01
## factor(vis_minority_2011.Q)3:pct_esp_vert_high_2011ct   0.20165 8.402e-01
## factor(vis_minority_2011.Q)4:pct_esp_vert_high_2011ct   3.09138 1.992e-03
## factor(vis_minority_2011.Q)5:pct_esp_vert_high_2011ct   4.57986 4.653e-06
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1378904 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.01063  
## # of obs: 688; # of groups: ct_no, 688 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -1797.177

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- clean_bei$df
.t["predict"] <- as.vector(predict(res.lmm))

g_vis_tc <- ggplot(.t, aes(y=100*predict/ct_area_hct, x=pct_esp_vert_high_2011ct, colour=factor(vis_minority_2011.Q), fill=factor(vis_minority_2011.Q))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='vis_minority_2011.Q', fill='vis_minority_2011.Q')
g_vis_tc

6.3 UI vs gentrified CT

No new models here, as gentrified variable is already used as a pure categorical variable.

7 Urban conditions and SES evolution

Tackling part 2 of objective #1: a multilevel model to test the reduction in socio-economic inequalities in urban conditions \[𝔼(UC_{ij} \mid X_{ij}) = \beta_{0j} + \beta_{1j} ∗ SES + \beta_{2j} ∗ Time + \beta_{3j} ∗ SES ∗ Time + \epsilon_{ij}\]

Data for bike lanes is available for 2006, 2011 and 2016 whereas greenness (and trees) is only available for 2011 and 2016.

# Get SCOREMAT
.bei_df_mlm.1 <- bei_df_aoi %>%
  select(CT_UID, starts_with("wSCOREMAT") & ends_with(".Q")) %>%
  pivot_longer(cols = starts_with("wSCOREMAT"), 
               names_to = "Year", names_pattern = "wSCOREMAT.(\\d{4}).Q",
               values_to = "wSCOREMAT.Q", values_drop_na = TRUE)
# Get vis_minority
.bei_df_mlm.2 <- bei_df_aoi %>%
  select(CT_UID, starts_with("vis_minority") & ends_with(".Q")) %>%
  pivot_longer(cols = starts_with("vis_minority"), 
               names_to = "Year", names_pattern = "vis_minority_(\\d{4}).Q",
               values_to = "vis_minority.Q", values_drop_na = TRUE)
# Get gentrified flag
.bei_df_mlm.3 <- bei_df_aoi %>%
  select(CT_UID, starts_with("gentrified_")) %>%
  pivot_longer(cols = starts_with("gentrified_"), 
               names_to = "Year_span", names_prefix = "gentrified_", 
               values_to = "gentrified", values_drop_na = TRUE) %>%
  extract(Year_span, "Year", "(\\d{4})_*")
# Get bike_lane
.bei_df_mlm.4 <- bei_df_aoi %>%
  select(CT_UID, starts_with("Bike_lane_total.hm")) %>%
  pivot_longer(cols = starts_with("Bike_lane_total.hm"), 
               names_to = "Year", names_pattern = "Bike_lane_total.hm.(\\d{4})ct",
               values_to = "Bike_lane_total.hm", values_drop_na = TRUE)
# Get Canopy
.bei_df_mlm.5 <- bei_df_aoi %>%
  select(CT_UID, starts_with("area_esp_vert_20")) %>%
  pivot_longer(cols = starts_with("area_esp_vert_20"), 
               names_to = "Year", names_pattern = "area_esp_vert_(\\d{4})ct",
               values_to = "area_esp_vert", values_drop_na = TRUE) %>%
  filter(Year %in% c('2011', '2017')) %>%
  mutate(Year = case_when(Year == '2017' ~ '2016',
                          TRUE ~ Year))
# Get Canopy
.bei_df_mlm.6 <- bei_df_aoi %>%
  select(CT_UID, starts_with("area_esp_vert_high_20")) %>%
  pivot_longer(cols = starts_with("area_esp_vert_high_20"), 
               names_to = "Year", names_pattern = "area_esp_vert_high_(\\d{4})ct",
               values_to = "area_esp_vert_high", values_drop_na = TRUE) %>%
  filter(Year %in% c('2011', '2017')) %>%
  mutate(Year = case_when(Year == '2017' ~ '2016',
                          TRUE ~ Year))

# Combine all subsets
bei_df_aoi_mlm <- bei_df_aoi %>%
  select(CT_UID, Population, street_length.hm, ct_area_hct) %>%
  full_join(.bei_df_mlm.1, by=c("CT_UID")) %>%
  full_join(.bei_df_mlm.2, by=c("CT_UID", "Year")) %>%
  full_join(.bei_df_mlm.3, by=c("CT_UID", "Year")) %>%
  full_join(.bei_df_mlm.4, by=c("CT_UID", "Year")) %>%
  full_join(.bei_df_mlm.5, by=c("CT_UID", "Year")) %>%
  full_join(.bei_df_mlm.6, by=c("CT_UID", "Year")) %>%
  units::drop_units() %>%
  #mutate(Year = factor(Year, levels = c("2006", "2011", "2016")))
  mutate(Yr = factor(Year, levels = c("2006", "2011", "2016")),
         Year = as.numeric(Year),
         Epoch = case_when(Yr == "2006" ~ 0,
                           Yr == "2011" ~ 1,
                           Yr == "2016" ~ 2))

# Get neighborhood matrix for the CT in MLM dataset
.CT16_mlm <- CT16 %>%
  transmute(CT_UID = GeoUID) %>%
  filter(CT_UID %in% unique(bei_df_aoi_mlm$CT_UID)) %>%
  mutate(ct_no = row_number())

# Compute neighborhood matrix
.ct_nb <- poly2nb(.CT16_mlm)
nbmatx <- nb2mat(.ct_nb, style = "B", zero.policy = TRUE)

# Add ct_no to mlm dataset
bei_df_aoi_mlm <- .CT16_mlm %>%
  as.data.frame() %>%
  select(CT_UID, ct_no) %>%
  inner_join(bei_df_aoi_mlm, by = "CT_UID")

7.1 UC vs Pampalon | material

7.1.1 Bike lane length

res.lmm <- fitme(Bike_lane_total.hm ~ wSCOREMAT.Q * Epoch + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm ~ wSCOREMAT.Q * Epoch + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                   Estimate Cond. SE  t-value   p-value
## (Intercept)       -2.80195 0.115135 -24.3361 0.000e+00
## wSCOREMAT.Q       -0.07789 0.015967  -4.8784 1.069e-06
## Epoch              0.24414 0.014855  16.4347 0.000e+00
## wSCOREMAT.Q:Epoch  0.00497 0.005228   0.9507 3.417e-01
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1339117 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  1.468  
## # of obs: 2072; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -6654.475

ci.lmm.1 <- confint(res.lmm, "wSCOREMAT.Q")

## lower wSCOREMAT.Q upper wSCOREMAT.Q 
##       -0.10921263       -0.04659658

ci.lmm.2 <- confint(res.lmm, "Epoch")

## lower Epoch upper Epoch 
##   0.2150454   0.2732732

ci.lmm.3 <- confint(res.lmm, "wSCOREMAT.Q:Epoch")

## lower wSCOREMAT.Q:Epoch upper wSCOREMAT.Q:Epoch 
##            -0.005272308             0.015218191

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm$data
.t["predict"] <- as.vector(predict(res.lmm))

ggplot(.t, aes(y=100*predict/street_length.hm, x=wSCOREMAT.Q, colour=factor(Year), fill=factor(Year))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='Year', fill='Year')

# AIC
AIC(res.lmm)

7.1.2 Bike lane length change (hurdle)

Bike lane length change

# Need to use PQL as ML and REML do not converge
res.lmm.0 <- fitme(I(Bike_lane_total.hm != 0) ~ wSCOREMAT.Q * Epoch + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = binomial(),
                 method = "PQL/L")
                 #verbose=c(TRACE=FALSE), control.HLfit = list(max.iter.mean = 500))
summary(res.lmm.0, details = c(p_value="Wald"))

## formula: I(Bike_lane_total.hm != 0) ~ wSCOREMAT.Q * Epoch + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by h-likelihood approximation.
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                   Estimate Cond. SE t-value   p-value
## (Intercept)       -3.54448  0.61826 -5.7330 9.867e-09
## wSCOREMAT.Q       -0.11669  0.15296 -0.7629 4.455e-01
## Epoch              1.95070  0.29917  6.5204 7.011e-11
## wSCOREMAT.Q:Epoch -0.03347  0.08496 -0.3939 6.936e-01
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1356573 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  13.68  
## # of obs: 2072; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                          logLik
##        h-likelihood: -1905.1069
## p_v(h) (marginal L):  -762.1575

# No canfidence for PQL method
# ci.lmm.0.1 <- confint(res.lmm.0, "wSCOREMAT.Q")
# ci.lmm.0.2 <- confint(res.lmm.0, "Epoch")
# ci.lmm.0.3 <- confint(res.lmm.0, "wSCOREMAT.Q:Epoch")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
  plot(simout)
}

# Plotting interaction || MEANINGLESS
# .t <- res.lmm.0$data
# .t["predict"] <- as.vector(predict(res.lmm.0))
# 
# ggplot(.t, aes(y=100*predict/street_length.hm, x=wSCOREMAT.Q, colour=factor(Year), fill=factor(Year))) + 
#   geom_smooth(method = "lm", alpha=.1) +
#   labs(colour='Year', fill='Year')

# Count model (Poisson)
bei_df_aoi_mlm_hurdle <- bei_df_aoi_mlm %>%
  select(-ct_no) %>%
  filter(Bike_lane_total.hm != 0)

# Get neighborhood matrix for the CT in MLM dataset
.CT16_mlm_hurdle <- CT16 %>%
  transmute(CT_UID = GeoUID) %>%
  filter(CT_UID %in% unique(bei_df_aoi_mlm_hurdle$CT_UID)) %>%
  mutate(ct_no = row_number())

# Compute neighborhood matrix
.ct_nb_hurdle <- poly2nb(.CT16_mlm_hurdle)
nbmatx_hurdle <- nb2mat(.ct_nb_hurdle, style = "B", zero.policy = TRUE)

# Add ct_no to mlm dataset
bei_df_aoi_mlm_hurdle <- .CT16_mlm_hurdle %>%
  as.data.frame() %>%
  select(CT_UID, ct_no) %>%
  inner_join(bei_df_aoi_mlm_hurdle, by = "CT_UID")


res.lmm.cnt <- fitme(Bike_lane_total.hm ~ wSCOREMAT.Q * Epoch + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm_hurdle, adjMatrix = nbmatx_hurdle,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm ~ wSCOREMAT.Q * Epoch + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                   Estimate Cond. SE t-value p-value
## (Intercept)       -2.30691 0.063544 -36.304 0.00000
## wSCOREMAT.Q       -0.03728 0.014976  -2.489 0.01279
## Epoch              0.22502 0.015110  14.892 0.00000
## wSCOREMAT.Q:Epoch -0.01592 0.005368  -2.966 0.00302
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1263186 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.4079  
## # of obs: 1432; # of groups: ct_no, 578 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -4878.605

# ci.lmm.cnt.1 <- confint(res.lmm.cnt, "wSCOREMAT.Q")
# ci.lmm.cnt.2 <- confint(res.lmm.cnt, "Epoch")
# ci.lmm.cnt.3 <- confint(res.lmm.cnt, "wSCOREMAT.Q:Epoch")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm.cnt$data
.t["predict"] <- as.vector(predict(res.lmm.cnt))

ggplot(.t, aes(y=100*predict/street_length.hm, x=wSCOREMAT.Q, colour=factor(Year), fill=factor(Year))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='Year', fill='Year')

# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)

# Combined AIC
"Combined Hurdle AIC:"

## [1] "Combined Hurdle AIC:"

-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)

## [1] 11293.52

7.1.3 Greenness

res.lmm <- fitme(area_esp_vert ~ wSCOREMAT.Q * Epoch + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert ~ wSCOREMAT.Q * Epoch + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                    Estimate Cond. SE t-value   p-value
## (Intercept)       -0.796231 0.037217 -21.394 0.000e+00
## wSCOREMAT.Q       -0.051170 0.010444  -4.899 9.608e-07
## Epoch              0.107861 0.012788   8.435 0.000e+00
## wSCOREMAT.Q:Epoch -0.006602 0.004752  -1.389 1.647e-01
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1360916 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.05332  
## # of obs: 1381; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -4569.472

ci.lmm.1 <- confint(res.lmm, "wSCOREMAT.Q")

## lower wSCOREMAT.Q upper wSCOREMAT.Q 
##       -0.07171847       -0.03056099

ci.lmm.2 <- confint(res.lmm, "Epoch")

## lower Epoch upper Epoch 
##  0.08280103  0.13292552

ci.lmm.3 <- confint(res.lmm, "wSCOREMAT.Q:Epoch")

## lower wSCOREMAT.Q:Epoch upper wSCOREMAT.Q:Epoch 
##            -0.015914020             0.002710266

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm$data
.t["predict"] <- as.vector(predict(res.lmm))

ggplot(.t, aes(y=100*predict/ct_area_hct, x=wSCOREMAT.Q, colour=factor(Year), fill=factor(Year))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='Year', fill='Year')

7.1.4 Canopy (trees)

res.lmm <- fitme(area_esp_vert_high ~ wSCOREMAT.Q * Epoch + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_high ~ wSCOREMAT.Q * Epoch + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                    Estimate Cond. SE t-value   p-value
## (Intercept)       -1.453075 0.057163 -25.420 0.000e+00
## wSCOREMAT.Q       -0.094496 0.016319  -5.790 7.022e-09
## Epoch              0.082641 0.019394   4.261 2.034e-05
## wSCOREMAT.Q:Epoch  0.008533 0.007427   1.149 2.506e-01
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1355721 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.1309  
## # of obs: 1381; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -3953.017

ci.lmm.1 <- confint(res.lmm, "wSCOREMAT.Q")

## lower wSCOREMAT.Q upper wSCOREMAT.Q 
##       -0.12663461       -0.06227973

ci.lmm.2 <- confint(res.lmm, "Epoch")

## lower Epoch upper Epoch 
##  0.04463734  0.12065176

ci.lmm.3 <- confint(res.lmm, "wSCOREMAT.Q:Epoch")

## lower wSCOREMAT.Q:Epoch upper wSCOREMAT.Q:Epoch 
##            -0.006016247             0.023086753

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm$data
.t["predict"] <- as.vector(predict(res.lmm))

ggplot(.t, aes(y=100*predict/ct_area_hct, x=wSCOREMAT.Q, colour=factor(Year), fill=factor(Year))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='Year', fill='Year')

7.2 UC vs Visible minority

7.2.1 Bike lane length

res.lmm <- fitme(Bike_lane_total.hm ~ vis_minority.Q * Epoch + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm ~ vis_minority.Q * Epoch + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                      Estimate Cond. SE t-value   p-value
## (Intercept)          -2.91941 0.121507 -24.027 0.000e+00
## vis_minority.Q       -0.04051 0.017908  -2.262 2.369e-02
## Epoch                 0.17212 0.020231   8.508 0.000e+00
## vis_minority.Q:Epoch  0.02563 0.005993   4.277 1.896e-05
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1338875 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  1.479  
## # of obs: 2065; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -6633.678

ci.lmm.1 <- confint(res.lmm, "vis_minority.Q")

## lower vis_minority.Q upper vis_minority.Q 
##         -0.075669803         -0.005351675

ci.lmm.2 <- confint(res.lmm, "Epoch")

## lower Epoch upper Epoch 
##   0.1324853   0.2117859

ci.lmm.3 <- confint(res.lmm, "vis_minority.Q:Epoch")

## lower vis_minority.Q:Epoch upper vis_minority.Q:Epoch 
##                 0.01388819                 0.03737936

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm$data
.t["predict"] <- as.vector(predict(res.lmm))

ggplot(.t, aes(y=100*predict/street_length.hm, x=vis_minority.Q, colour=factor(Year), fill=factor(Year))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='Year', fill='Year')

# AIC
AIC(res.lmm)

7.2.2 Bike lane length change (hurdle)

Bike lane length change

# Need to use PQL as ML and REML do not converge
res.lmm.0 <- fitme(I(Bike_lane_total.hm != 0) ~ vis_minority.Q * Epoch + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = binomial(),
                 method = "PQL/L")
                 #verbose=c(TRACE=FALSE), control.HLfit = list(max.iter.mean = 500))
summary(res.lmm.0, details = c(p_value="Wald"))

## formula: I(Bike_lane_total.hm != 0) ~ vis_minority.Q * Epoch + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by h-likelihood approximation.
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                      Estimate Cond. SE t-value   p-value
## (Intercept)          -2.82256   0.8050 -3.5063 4.544e-04
## vis_minority.Q       -0.32008   0.1995 -1.6046 1.086e-01
## Epoch                 2.13641   0.4621  4.6229 3.785e-06
## vis_minority.Q:Epoch -0.07509   0.1182 -0.6355 5.251e-01
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1355477 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  13.89  
## # of obs: 2065; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                          logLik
##        h-likelihood: -1905.7051
## p_v(h) (marginal L):  -756.7773

# No canfidence for PQL method
# ci.lmm.0.1 <- confint(res.lmm.0, "wSCOREMAT.Q")
# ci.lmm.0.2 <- confint(res.lmm.0, "Epoch")
# ci.lmm.0.3 <- confint(res.lmm.0, "wSCOREMAT.Q:Epoch")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
  plot(simout)
}

# Plotting interaction || MEANINGLESS
# .t <- res.lmm.0$data
# .t["predict"] <- as.vector(predict(res.lmm.0))
# 
# ggplot(.t, aes(y=100*predict/street_length.hm, x=vis_minority.Q, colour=factor(Year), fill=factor(Year))) + 
#   geom_smooth(method = "lm", alpha=.1) +
#   labs(colour='Year', fill='Year')

# Count model (Poisson)
bei_df_aoi_mlm_hurdle <- bei_df_aoi_mlm %>%
  select(-ct_no) %>%
  filter(Bike_lane_total.hm != 0)

# Get neighborhood matrix for the CT in MLM dataset
.CT16_mlm_hurdle <- CT16 %>%
  transmute(CT_UID = GeoUID) %>%
  filter(CT_UID %in% unique(bei_df_aoi_mlm_hurdle$CT_UID)) %>%
  mutate(ct_no = row_number())

# Compute neighborhood matrix
.ct_nb_hurdle <- poly2nb(.CT16_mlm_hurdle)
nbmatx_hurdle <- nb2mat(.ct_nb_hurdle, style = "B", zero.policy = TRUE)

# Add ct_no to mlm dataset
bei_df_aoi_mlm_hurdle <- .CT16_mlm_hurdle %>%
  as.data.frame() %>%
  select(CT_UID, ct_no) %>%
  inner_join(bei_df_aoi_mlm_hurdle, by = "CT_UID")


res.lmm.cnt <- fitme(Bike_lane_total.hm ~ vis_minority.Q * Epoch + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm_hurdle, adjMatrix = nbmatx_hurdle,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm ~ vis_minority.Q * Epoch + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                       Estimate Cond. SE  t-value p-value
## (Intercept)          -2.471864 0.075300 -32.8268  0.0000
## vis_minority.Q        0.011495 0.017116   0.6716  0.5018
## Epoch                 0.196212 0.020512   9.5658  0.0000
## vis_minority.Q:Epoch -0.004416 0.006126  -0.7209  0.4710
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##    1.rho 
## 0.128713 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.4092  
## # of obs: 1426; # of groups: ct_no, 578 
##  ------------- Likelihood values  -------------
##                        logLik
## p_v(h) (marginal L): -4869.73

# ci.lmm.cnt.1 <- confint(res.lmm.cnt, "wSCOREMAT.Q")
# ci.lmm.cnt.2 <- confint(res.lmm.cnt, "Epoch")
# ci.lmm.cnt.3 <- confint(res.lmm.cnt, "wSCOREMAT.Q:Epoch")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm.cnt$data
.t["predict"] <- as.vector(predict(res.lmm.cnt))

ggplot(.t, aes(y=100*predict/street_length.hm, x=vis_minority.Q, colour=factor(Year), fill=factor(Year))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='Year', fill='Year')

# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)

# Combined AIC
"Combined Hurdle AIC:"

## [1] "Combined Hurdle AIC:"

-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)

## [1] 11265.01

7.2.3 Greenness

res.lmm <- fitme(area_esp_vert ~ vis_minority.Q * Epoch + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert ~ vis_minority.Q * Epoch + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                      Estimate Cond. SE t-value   p-value
## (Intercept)          -0.90488 0.046273 -19.555 0.000e+00
## vis_minority.Q       -0.01317 0.012300  -1.071 2.843e-01
## Epoch                 0.18691 0.015869  11.778 0.000e+00
## vis_minority.Q:Epoch -0.03039 0.004826  -6.297 3.035e-10
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##    1.rho 
## 0.136075 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.0568  
## # of obs: 1378; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -4557.123

ci.lmm.1 <- confint(res.lmm, "vis_minority.Q")

## lower vis_minority.Q upper vis_minority.Q 
##          -0.03735459           0.01112046

ci.lmm.2 <- confint(res.lmm, "Epoch")

## lower Epoch upper Epoch 
##   0.1558138   0.2180168

ci.lmm.3 <- confint(res.lmm, "vis_minority.Q:Epoch")

## lower vis_minority.Q:Epoch upper vis_minority.Q:Epoch 
##                -0.03984867                -0.02093196

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm$data
.t["predict"] <- as.vector(predict(res.lmm))

ggplot(.t, aes(y=100*predict/ct_area_hct, x=vis_minority.Q, colour=factor(Year), fill=factor(Year))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='Year', fill='Year')

7.2.4 Canopy (trees)

res.lmm <- fitme(area_esp_vert_high ~ vis_minority.Q * Epoch + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_high ~ vis_minority.Q * Epoch + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                       Estimate Cond. SE t-value   p-value
## (Intercept)          -1.390055 0.070353 -19.758 0.000e+00
## vis_minority.Q       -0.106805 0.018973  -5.629 1.809e-08
## Epoch                 0.081427 0.024960   3.262 1.105e-03
## vis_minority.Q:Epoch  0.009049 0.007679   1.178 2.387e-01
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1355641 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.1351  
## # of obs: 1378; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                        logLik
## p_v(h) (marginal L): -3949.38

ci.lmm.1 <- confint(res.lmm, "vis_minority.Q")

## lower vis_minority.Q upper vis_minority.Q 
##          -0.14411501          -0.06937289

ci.lmm.2 <- confint(res.lmm, "Epoch")

## lower Epoch upper Epoch 
##  0.03251626  0.13034790

ci.lmm.3 <- confint(res.lmm, "vis_minority.Q:Epoch")

## lower vis_minority.Q:Epoch upper vis_minority.Q:Epoch 
##               -0.005999146                0.024098267

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm$data
.t["predict"] <- as.vector(predict(res.lmm))

ggplot(.t, aes(y=100*predict/ct_area_hct, x=vis_minority.Q, colour=factor(Year), fill=factor(Year))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='Year', fill='Year')

7.3 UC vs gentrified CT

7.3.1 Bike lane length

res.lmm <- fitme(Bike_lane_total.hm ~ gentrified * Epoch + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm ~ gentrified * Epoch + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                      Estimate Cond. SE t-value   p-value
## (Intercept)           -2.9821 0.105790 -28.189 0.000e+00
## gentrifiedTRUE        -0.2936 0.037448  -7.839 4.552e-15
## Epoch                  0.2169 0.008302  26.128 0.000e+00
## gentrifiedTRUE:Epoch   0.1918 0.023635   8.115 4.441e-16
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1339985 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  1.464  
## # of obs: 1997; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -6344.132

ci.lmm.1 <- confint(res.lmm, "gentrifiedTRUE")

## lower gentrifiedTRUE upper gentrifiedTRUE 
##           -0.3672111           -0.2204272

ci.lmm.2 <- confint(res.lmm, "Epoch")

## lower Epoch upper Epoch 
##   0.2006498   0.2331913

ci.lmm.3 <- confint(res.lmm, "gentrifiedTRUE:Epoch")

## lower gentrifiedTRUE:Epoch upper gentrifiedTRUE:Epoch 
##                  0.1455700                  0.2382088

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm$data
.t["predict"] <- as.vector(predict(res.lmm))

#ggplot(.t, aes(y=100*predict/street_length.hm, x=gentrified, colour=factor(Year), fill=factor(Year))) + 
#  geom_boxplot(alpha = .1) +
#  labs(colour='Year', fill='Year')

ggplot(.t, aes(y=100*predict/street_length.hm, x=factor(Year), colour=factor(gentrified), fill=factor(gentrified))) + 
  geom_boxplot(alpha = .1) +
  labs(colour='gentrified', fill='gentrified')

# AIC
AIC(res.lmm)

7.3.2 Bike lane length change (hurdle)

Bike lane length change

# Need to use PQL as ML and REML do not converge
res.lmm.0 <- fitme(I(Bike_lane_total.hm != 0) ~ gentrified * Epoch + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = binomial(),
                 method = "PQL/L")
                 #verbose=c(TRACE=FALSE), control.HLfit = list(max.iter.mean = 500))
summary(res.lmm.0, details = c(p_value="Wald"))

## formula: I(Bike_lane_total.hm != 0) ~ gentrified * Epoch + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by h-likelihood approximation.
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: binomial( link = logit ) 
##  ------------ Fixed effects (beta) ------------
##                      Estimate Cond. SE t-value p-value
## (Intercept)           -3.6848   0.4257  -8.656 0.00000
## gentrifiedTRUE        -0.8627   0.3831  -2.252 0.02432
## Epoch                  1.7501   0.1722  10.164 0.00000
## gentrifiedTRUE:Epoch   0.4292   0.2957   1.451 0.14666
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1356131 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  14.86  
## # of obs: 1997; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                          logLik
##        h-likelihood: -1908.9928
## p_v(h) (marginal L):  -732.2444

# No canfidence for PQL method
# ci.lmm.0.1 <- confint(res.lmm.0, "wSCOREMAT.Q")
# ci.lmm.0.2 <- confint(res.lmm.0, "Epoch")
# ci.lmm.0.3 <- confint(res.lmm.0, "wSCOREMAT.Q:Epoch")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.0, plot = F)
  plot(simout)
}

# Plotting interaction || MEANINGLESS
# .t <- res.lmm.0$data
# .t["predict"] <- as.vector(predict(res.lmm.0))
# 
# ggplot(.t, aes(y=100*predict/street_length.hm, x=gentrified, colour=factor(Year), fill=factor(Year))) + 
#   geom_smooth(method = "lm", alpha=.1) +
#   labs(colour='Year', fill='Year')

# Count model (Poisson)
bei_df_aoi_mlm_hurdle <- bei_df_aoi_mlm %>%
  select(-ct_no) %>%
  filter(Bike_lane_total.hm != 0)

# Get neighborhood matrix for the CT in MLM dataset
.CT16_mlm_hurdle <- CT16 %>%
  transmute(CT_UID = GeoUID) %>%
  filter(CT_UID %in% unique(bei_df_aoi_mlm_hurdle$CT_UID)) %>%
  mutate(ct_no = row_number())

# Compute neighborhood matrix
.ct_nb_hurdle <- poly2nb(.CT16_mlm_hurdle)
nbmatx_hurdle <- nb2mat(.ct_nb_hurdle, style = "B", zero.policy = TRUE)

# Add ct_no to mlm dataset
bei_df_aoi_mlm_hurdle <- .CT16_mlm_hurdle %>%
  as.data.frame() %>%
  select(CT_UID, ct_no) %>%
  inner_join(bei_df_aoi_mlm_hurdle, by = "CT_UID")


res.lmm.cnt <- fitme(Bike_lane_total.hm ~ gentrified * Epoch + offset(log(street_length.hm)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm_hurdle, adjMatrix = nbmatx_hurdle,
                 family = Tpoisson())
summary(res.lmm.cnt, details = c(p_value="Wald"))

## formula: Bike_lane_total.hm ~ gentrified * Epoch + offset(log(street_length.hm)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: 0-truncated poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                      Estimate Cond. SE t-value p-value
## (Intercept)          -2.41337 0.048693 -49.563 0.00000
## gentrifiedTRUE       -0.05786 0.038841  -1.490 0.13631
## Epoch                 0.16722 0.008414  19.874 0.00000
## gentrifiedTRUE:Epoch  0.04232 0.024442   1.731 0.08338
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1292178 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.4032  
## # of obs: 1377; # of groups: ct_no, 578 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -4672.713

# ci.lmm.cnt.1 <- confint(res.lmm.cnt, "wSCOREMAT.Q")
# ci.lmm.cnt.2 <- confint(res.lmm.cnt, "Epoch")
# ci.lmm.cnt.3 <- confint(res.lmm.cnt, "wSCOREMAT.Q:Epoch")

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm.cnt, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm.cnt$data
.t["predict"] <- as.vector(predict(res.lmm.cnt))

ggplot(.t, aes(y=100*predict/street_length.hm, x=gentrified, colour=factor(Year), fill=factor(Year))) + 
  geom_smooth(method = "lm", alpha=.1) +
  labs(colour='Year', fill='Year')

# AIC
AIC(res.lmm.0)
AIC(res.lmm.cnt)

# Combined AIC
"Combined Hurdle AIC:"

## [1] "Combined Hurdle AIC:"

-2 * res.lmm.0$APHLs$p_v - 2 * res.lmm.cnt$APHLs$p_v + 2 * Reduce(sum, res.lmm.cnt$dfs)

## [1] 10821.91

7.3.3 Greenness

res.lmm <- fitme(area_esp_vert ~ gentrified * Epoch + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert ~ gentrified * Epoch + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                      Estimate Cond. SE  t-value p-value
## (Intercept)          -0.94510  0.02555 -36.9848 0.00000
## gentrifiedTRUE       -0.07346  0.04196  -1.7506 0.08001
## Epoch                 0.08634  0.00650  13.2828 0.00000
## gentrifiedTRUE:Epoch  0.01826  0.02673   0.6831 0.49452
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1360636 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.05927  
## # of obs: 1363; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -4548.274

ci.lmm.1 <- confint(res.lmm, "gentrifiedTRUE")

## lower gentrifiedTRUE upper gentrifiedTRUE 
##         -0.155758976          0.008712541

ci.lmm.2 <- confint(res.lmm, "Epoch")

## lower Epoch upper Epoch 
##  0.07360361  0.09908424

ci.lmm.3 <- confint(res.lmm, "gentrifiedTRUE:Epoch")

## lower gentrifiedTRUE:Epoch upper gentrifiedTRUE:Epoch 
##                -0.03416116                 0.07070056

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm$data
.t["predict"] <- as.vector(predict(res.lmm))

#ggplot(.t, aes(y=100*predict/street_length.hm, x=gentrified, colour=factor(Year), fill=factor(Year))) + 
#  geom_boxplot(alpha = .1) +
#  labs(colour='Year', fill='Year')

ggplot(.t, aes(y=100*predict/ct_area_hct, x=factor(Year), colour=factor(gentrified), fill=factor(gentrified))) + 
  geom_boxplot(alpha = .1) +
  labs(colour='gentrified', fill='gentrified')

7.3.4 Canopy (trees)

res.lmm <- fitme(area_esp_vert_high ~ gentrified * Epoch + offset(log(ct_area_hct)) + adjacency(1|ct_no),
                 data = bei_df_aoi_mlm, adjMatrix = nbmatx,
                 family = poisson())
summary(res.lmm, details = c(p_value="Wald"))

## formula: area_esp_vert_high ~ gentrified * Epoch + offset(log(ct_area_hct)) + 
##     adjacency(1 | ct_no)
## Estimation of corrPars and lambda by Laplace ML approximation (p_v).
## Estimation of fixed effects by Laplace ML approximation (p_v).
## Estimation of lambda by 'outer' ML, maximizing p_v.
## family: poisson( link = log ) 
##  ------------ Fixed effects (beta) ------------
##                      Estimate Cond. SE t-value  p-value
## (Intercept)          -1.71058  0.03924 -43.593 0.000000
## gentrifiedTRUE       -0.18037  0.06354  -2.839 0.004531
## Epoch                 0.08906  0.01008   8.834 0.000000
## gentrifiedTRUE:Epoch  0.09291  0.03944   2.355 0.018500
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##     1.rho 
## 0.1355664 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    ct_no  :  0.1443  
## # of obs: 1363; # of groups: ct_no, 705 
##  ------------- Likelihood values  -------------
##                         logLik
## p_v(h) (marginal L): -3930.354

ci.lmm.1 <- confint(res.lmm, "gentrifiedTRUE")

## lower gentrifiedTRUE upper gentrifiedTRUE 
##          -0.30508466          -0.05604759

ci.lmm.2 <- confint(res.lmm, "Epoch")

## lower Epoch upper Epoch 
##   0.0693051   0.1088227

ci.lmm.3 <- confint(res.lmm, "gentrifiedTRUE:Epoch")

## lower gentrifiedTRUE:Epoch upper gentrifiedTRUE:Epoch 
##                  0.0156415                  0.1702350

if (plot.resid) {
  simout <- simulateResiduals(fittedModel = res.lmm, plot = F)
  plot(simout)
}

# Plotting interaction
.t <- res.lmm$data
.t["predict"] <- as.vector(predict(res.lmm))

ggplot(.t, aes(y=100*predict/street_length.hm, x=gentrified, colour=factor(Year), fill=factor(Year))) + 
  geom_boxplot(alpha = .1) +
  labs(colour='Year', fill='Year')

ggplot(.t, aes(y=100*predict/ct_area_hct, x=factor(Year), colour=factor(gentrified), fill=factor(gentrified))) + 
  geom_boxplot(alpha = .1) +
  labs(colour='gentrified', fill='gentrified')

8 Summary of association between SES and Built Environment

NOTICE: This section needs to be updated

Presenting interaction terms (STROBE and more specifically here): best practice suggests we should include interaction terms for specific values of covariates.

SES metric	BE metric	UC 2011 trend	UI 2011 -> 2016 trend (SES effect only, no interaction)
Pampalon / MAT 2011	Bike lane	0.857 (0.781 – 0.94)	0.934 (0.882 – 0.988)
Pampalon / MAT 2011	Greenness	0.92 (0.902 – 0.939)	0.944 (0.918 – 0.971)
Pampalon / MAT 2011	Tree canopy	0.887 (0.86 – 0.914)	0.892 (0.865 – 0.919)
Visible Minority 2011	Bike lane	0.856 (0.759 – 0.965)	0.894 (0.829 – 0.963)
Visible Minority 2011	Greenness	0.914 (0.89 – 0.938)	0.918 (0.887 – 0.95)
Visible Minority 2011	Tree canopy	0.87 (0.836 – 0.906)	0.903 (0.87 – 0.937)
Gentrified 2011-2016	Bike lane	0.936 (0.717 – 1.22)	1.08 (0.905 – 1.286)
Gentrified 2011-2016	Greenness	0.822 (0.769 – 0.878)	0.817 (0.719 – 0.927)
Gentrified 2011-2016	Tree canopy	0.86 (0.78 – 0.949)	0.815 (0.708 – 0.937)

Italic: not significant

SES metric	BE metric	UC 2011 trend	UI 2011 -> 2016 trend (SES effect only, no interaction)
Pampalon / MAT 2011 (hurdle + cnt)	Bike lane	0.810 (p=0.016) & 0.941 (p=0.021)	0.98 (p=0.9205) & 0.984 (p=0.4830)
Pampalon / MAT 2011	Greenness	0.92 (p<0.001)	0.944 (p<0.001)
Pampalon / MAT 2011	Tree canopy	0.887 (p<0.001)	0.892 (p<0.001)
Visible Minority 2011 (hurdle + cnt)	Bike lane	0.801 (p=0.054*) & 0.897 (p=0.001)	0.59 (p=0.0771) & 0.967 (p=0.259)
Visible Minority 2011	Greenness	0.914 (p<0.001)	0.918 (p<0.001)
Visible Minority 2011	Tree canopy	0.87 (p<0.001)	0.903 (p<0.001)
Gentrified 2011-2016 (hurdle + cnt)	Bike lane	0.936 (p=0.789) & 1.1 (p=0.236)	1.211 (p=0.747) & 1.224 (p=0.005)
Gentrified 2011-2016	Greenness	0.822 (p<0.001)	0.817 (p=0.002)
Gentrified 2011-2016	Tree canopy	0.86 (p=0.003)	0.815 (p=0.004)

p-values: Wald t-tests except when p-value close to threshold where Likelihood ratio test has been used (marked with *) Bold: significant

9 Tables and figures

9.1 SES and BEI outcomes

# Extract SES and BEI in 2011
.tab1_2011 <- bei_df_aoi %>%
  transmute(year = "2011",
            zone = zone,
            CT_UID = CT_UID,
            Population = Population,
            `Material score` = wSCOREMAT.2011,
            `% visible minorities` = vis_minority_2011,
            `gentrified` = gentrified_2016_2011,
            `Bike lane` = Bike_lane_total.2011ct > 0,
            `Bike lane length (hm)` = Bike_lane_total.2011ct/100,
            `% Bike lane to streets` = Bike_lane.by.street.2011ct,
            `Green space area (ha)` = area_esp_vert_2011ct,
            `% green space within CT` = pct_esp_vert_2011ct,
            `Tree canopy area (ha)` = area_esp_vert_high_2011ct,
            `% tree canopy within CT` = pct_esp_vert_high_2011ct)

# Extract SES and BEI in 2016
.tab1_2016 <- bei_df_aoi %>%
  transmute(year = "2016",
            zone = zone,
            CT_UID = CT_UID,
            `Bike lane` = Bike_lane_total.2016ct > 0,
            `Bike lane length (hm)` = Bike_lane_total.2016ct/100,
            `% Bike lane to streets` = Bike_lane.by.street.2016ct,
            `Green space area (ha)` = area_esp_vert_2017ct,
            `% green space within CT` = pct_esp_vert_2017ct,
            `Tree canopy area (ha)` = area_esp_vert_high_2017ct,
            `% tree canopy within CT` = pct_esp_vert_high_2017ct)

.tab1_med <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year) %>%
  summarize(across(where(is.numeric), median, na.rm=TRUE), 
            `% gentrified` = 100 * sum(`gentrified`, na.rm = TRUE) / n(),
            `Bike lane presence` = sum(`Bike lane`))
.tab1_p25 <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year) %>%
  summarize(across(where(is.numeric), quantile, probs=.25, na.rm=TRUE))
.tab1_p75 <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year) %>%
  summarize(across(where(is.numeric), quantile, probs=.75, na.rm=TRUE))
tab1 <- inner_join(.tab1_med, .tab1_p25, by="year", suffix = c("", " p25")) %>%
  inner_join(.tab1_p75, by="year", suffix = c("", " p75")) %>%
  transmute(year = year,
            Population = case_when(year == "2016" ~ "", 
                                   TRUE ~ paste0(round(`Population`,1), " (", round(`Population p25`,1), " - ", round(`Population p75`,1), ")")),
            `Material score` = case_when(year == "2016" ~ "", 
                                         TRUE ~ paste0(round(`Material score`,3), " (", round(`Material score p25`,3), " - ", round(`Material score p75`,3), ")")),
            `% visible minorities` = case_when(is.na(`% visible minorities`) ~ "", 
                                               TRUE ~ paste0(round(`% visible minorities`,1), " (", round(`% visible minorities p25`,1), " - ", round(`% visible minorities p75`,1), ")")),
            `% gentrified CTs` = case_when(year == "2016" ~ "", 
                                           TRUE ~ as.character(round(`% gentrified`, 1))),
            `Bike lane presence (N)` = as.character(`Bike lane presence`),
            `Bike lane length (hm)` = paste0(round(`Bike lane length (hm)`,1), " (", round(`Bike lane length (hm) p25`,1), " - ", round(`Bike lane length (hm) p75`,1), ")"),
            `% Bike lane to streets` = paste0(round(`% Bike lane to streets`,1), " (", round(`% Bike lane to streets p25`,1), " - ", round(`% Bike lane to streets p75`,1), ")"),
            `Green space area (ha)` = paste0(round(`Green space area (ha)`,1), " (", round(`Green space area (ha) p25`,1), " - ", round(`Green space area (ha) p75`,1), ")"),
            `% green space within CT` = paste0(round(`% green space within CT`,1), " (", round(`% green space within CT p25`,1), " - ", round(`% green space within CT p75`,1), ")"),
            `Tree canopy area (ha)` = paste0(round(`Tree canopy area (ha)`,1), " (", round(`Tree canopy area (ha) p25`,1), " - ", round(`Tree canopy area (ha) p75`,1), ")"),
            `% tree canopy within CT` = paste0(round(`% tree canopy within CT`,1), " (", round(`% tree canopy within CT p25`,1), " - ", round(`% tree canopy within CT p75`,1), ")")) %>%
  tibble::column_to_rownames("year") %>%
  t() %>% 
  as.data.frame()

tab1

9.1.1 By Q1 / Q5

We look at Urban conditions in 2011 and 2016 within Q1 and Q5 of 2011 equity metrics

# Extract SES and BEI in 2011
.tab1_2011 <- bei_df_aoi %>%
  filter(wSCOREMAT.2011.Q %in% c(1, 5) ) %>%
  transmute(year = "2011",
            quintile = case_when(wSCOREMAT.2011.Q == 1 ~ "Q1 Deprivation",
                                 wSCOREMAT.2011.Q == 5 ~ "Q5 Deprivation"),
            zone = zone,
            CT_UID = CT_UID,
            Population = Population,
            `Material score` = wSCOREMAT.2011,
            `% visible minorities` = vis_minority_2011,
            `gentrified` = gentrified_2016_2011,
            `Bike lane` = Bike_lane_total.2011ct > 0,
            `Bike lane length (hm)` = Bike_lane_total.2011ct/100,
            `% Bike lane to streets` = Bike_lane.by.street.2011ct,
            `Green space area (ha)` = area_esp_vert_2011ct,
            `% green space within CT` = pct_esp_vert_2011ct,
            `Tree canopy area (ha)` = area_esp_vert_high_2011ct,
            `% tree canopy within CT` = pct_esp_vert_high_2011ct)

# Extract SES and BEI in 2016
.tab1_2016 <- bei_df_aoi %>%
  filter(wSCOREMAT.2011.Q %in% c(1, 5) ) %>%
  transmute(year = "2016",
            quintile = case_when(wSCOREMAT.2011.Q == 1 ~ "Q1 Deprivation",
                                 wSCOREMAT.2011.Q == 5 ~ "Q5 Deprivation"),
            zone = zone,
            CT_UID = CT_UID,
            `Bike lane` = Bike_lane_total.2016ct > 0,
            `Bike lane length (hm)` = Bike_lane_total.2016ct/100,
            `% Bike lane to streets` = Bike_lane.by.street.2016ct,
            `Green space area (ha)` = area_esp_vert_2017ct,
            `% green space within CT` = pct_esp_vert_2017ct,
            `Tree canopy area (ha)` = area_esp_vert_high_2017ct,
            `% tree canopy within CT` = pct_esp_vert_high_2017ct)

.tab1_med <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year, quintile) %>%
  summarize(across(where(is.numeric), median, na.rm=TRUE), 
            `% gentrified` = 100 * sum(`gentrified`, na.rm = TRUE) / n(),
            `Bike lane presence` = sum(`Bike lane`))

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

.tab1_p25 <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year, quintile) %>%
  summarize(across(where(is.numeric), quantile, probs=.25, na.rm=TRUE))

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

.tab1_p75 <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year, quintile) %>%
  summarize(across(where(is.numeric), quantile, probs=.75, na.rm=TRUE))

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

tab1 <- inner_join(.tab1_med, .tab1_p25, by=c("year", "quintile"), suffix = c("", " p25")) %>%
  inner_join(.tab1_p75, by=c("year", "quintile"), suffix = c("", " p75")) %>%
  transmute(year.quintile = paste(year, quintile),
            Population = case_when(year == "2016" ~ "", 
                                   TRUE ~ paste0(round(`Population`,1), " (", round(`Population p25`,1), " - ", round(`Population p75`,1), ")")),
            `Material score` = case_when(year == "2016" ~ "", 
                                         TRUE ~ paste0(round(`Material score`,3), " (", round(`Material score p25`,3), " - ", round(`Material score p75`,3), ")")),
            `% visible minorities` = case_when(is.na(`% visible minorities`) ~ "", 
                                               TRUE ~ paste0(round(`% visible minorities`,1), " (", round(`% visible minorities p25`,1), " - ", round(`% visible minorities p75`,1), ")")),
            `% gentrified CTs` = case_when(year == "2016" ~ "", 
                                           TRUE ~ as.character(round(`% gentrified`, 1))),
            `Bike lane presence (N)` = as.character(`Bike lane presence`),
            `Bike lane length (hm)` = paste0(round(`Bike lane length (hm)`,1), " (", round(`Bike lane length (hm) p25`,1), " - ", round(`Bike lane length (hm) p75`,1), ")"),
            `% Bike lane to streets` = paste0(round(`% Bike lane to streets`,1), " (", round(`% Bike lane to streets p25`,1), " - ", round(`% Bike lane to streets p75`,1), ")"),
            `Green space area (ha)` = paste0(round(`Green space area (ha)`,1), " (", round(`Green space area (ha) p25`,1), " - ", round(`Green space area (ha) p75`,1), ")"),
            `% green space within CT` = paste0(round(`% green space within CT`,1), " (", round(`% green space within CT p25`,1), " - ", round(`% green space within CT p75`,1), ")"),
            `Tree canopy area (ha)` = paste0(round(`Tree canopy area (ha)`,1), " (", round(`Tree canopy area (ha) p25`,1), " - ", round(`Tree canopy area (ha) p75`,1), ")"),
            `% tree canopy within CT` = paste0(round(`% tree canopy within CT`,1), " (", round(`% tree canopy within CT p25`,1), " - ", round(`% tree canopy within CT p75`,1), ")")) %>%
  tibble::column_to_rownames("year.quintile") %>%
  t() %>% 
  as.data.frame()

tab1

# Extract SES and BEI in 2011
.tab1_2011 <- bei_df_aoi %>%
  filter(vis_minority_2011.Q %in% c(1, 5) ) %>%
  transmute(year = "2011",
            quintile = case_when(vis_minority_2011.Q == 1 ~ "Q1 Vis. Minorities",
                                 vis_minority_2011.Q == 5 ~ "Q5 Vis. Minorities"),
            zone = zone,
            CT_UID = CT_UID,
            Population = Population,
            `Material score` = wSCOREMAT.2011,
            `% visible minorities` = vis_minority_2011,
            `gentrified` = gentrified_2016_2011,
            `Bike lane` = Bike_lane_total.2011ct > 0,
            `Bike lane length (hm)` = Bike_lane_total.2011ct/100,
            `% Bike lane to streets` = Bike_lane.by.street.2011ct,
            `Green space area (ha)` = area_esp_vert_2011ct,
            `% green space within CT` = pct_esp_vert_2011ct,
            `Tree canopy area (ha)` = area_esp_vert_high_2011ct,
            `% tree canopy within CT` = pct_esp_vert_high_2011ct)

# Extract SES and BEI in 2016
.tab1_2016 <- bei_df_aoi %>%
  filter(vis_minority_2011.Q %in% c(1, 5) ) %>%
  transmute(year = "2016",
            quintile = case_when(vis_minority_2011.Q == 1 ~ "Q1 Vis. Minorities",
                                 vis_minority_2011.Q == 5 ~ "Q5 Vis. Minorities"),
            zone = zone,
            CT_UID = CT_UID,
            `Bike lane` = Bike_lane_total.2016ct > 0,
            `Bike lane length (hm)` = Bike_lane_total.2016ct/100,
            `% Bike lane to streets` = Bike_lane.by.street.2016ct,
            `Green space area (ha)` = area_esp_vert_2017ct,
            `% green space within CT` = pct_esp_vert_2017ct,
            `Tree canopy area (ha)` = area_esp_vert_high_2017ct,
            `% tree canopy within CT` = pct_esp_vert_high_2017ct)

.tab1_med <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year, quintile) %>%
  summarize(across(where(is.numeric), median, na.rm=TRUE), 
            `% gentrified` = 100 * sum(`gentrified`, na.rm = TRUE) / n(),
            `Bike lane presence` = sum(`Bike lane`))

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

.tab1_p25 <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year, quintile) %>%
  summarize(across(where(is.numeric), quantile, probs=.25, na.rm=TRUE))

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

.tab1_p75 <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year, quintile) %>%
  summarize(across(where(is.numeric), quantile, probs=.75, na.rm=TRUE))

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

tab1 <- inner_join(.tab1_med, .tab1_p25, by=c("year", "quintile"), suffix = c("", " p25")) %>%
  inner_join(.tab1_p75, by=c("year", "quintile"), suffix = c("", " p75")) %>%
  transmute(year.quintile = paste(year, quintile),
            Population = case_when(year == "2016" ~ "", 
                                   TRUE ~ paste0(round(`Population`,1), " (", round(`Population p25`,1), " - ", round(`Population p75`,1), ")")),
            `Material score` = case_when(year == "2016" ~ "", 
                                         TRUE ~ paste0(round(`Material score`,3), " (", round(`Material score p25`,3), " - ", round(`Material score p75`,3), ")")),
            `% visible minorities` = case_when(is.na(`% visible minorities`) ~ "", 
                                               TRUE ~ paste0(round(`% visible minorities`,1), " (", round(`% visible minorities p25`,1), " - ", round(`% visible minorities p75`,1), ")")),
            `% gentrified CTs` = case_when(year == "2016" ~ "", 
                                           TRUE ~ as.character(round(`% gentrified`, 1))),
            `Bike lane presence (N)` = as.character(`Bike lane presence`),
            `Bike lane length (hm)` = paste0(round(`Bike lane length (hm)`,1), " (", round(`Bike lane length (hm) p25`,1), " - ", round(`Bike lane length (hm) p75`,1), ")"),
            `% Bike lane to streets` = paste0(round(`% Bike lane to streets`,1), " (", round(`% Bike lane to streets p25`,1), " - ", round(`% Bike lane to streets p75`,1), ")"),
            `Green space area (ha)` = paste0(round(`Green space area (ha)`,1), " (", round(`Green space area (ha) p25`,1), " - ", round(`Green space area (ha) p75`,1), ")"),
            `% green space within CT` = paste0(round(`% green space within CT`,1), " (", round(`% green space within CT p25`,1), " - ", round(`% green space within CT p75`,1), ")"),
            `Tree canopy area (ha)` = paste0(round(`Tree canopy area (ha)`,1), " (", round(`Tree canopy area (ha) p25`,1), " - ", round(`Tree canopy area (ha) p75`,1), ")"),
            `% tree canopy within CT` = paste0(round(`% tree canopy within CT`,1), " (", round(`% tree canopy within CT p25`,1), " - ", round(`% tree canopy within CT p75`,1), ")")) %>%
  tibble::column_to_rownames("year.quintile") %>%
  t() %>% 
  as.data.frame()

tab1

9.1.2 By gentrification status

# Extract SES and BEI in 2011
.tab1_2011 <- bei_df_aoi %>%
  transmute(year = "2011",
            quintile = case_when(gentrified_2016_2011 ~ "Gentrified",
                                 TRUE ~ "Not gentrified"),
            zone = zone,
            CT_UID = CT_UID,
            Population = Population,
            `Material score` = wSCOREMAT.2011,
            `% visible minorities` = vis_minority_2011,
            `gentrified` = gentrified_2016_2011,
            `Bike lane` = Bike_lane_total.2011ct > 0,
            `Bike lane length (hm)` = Bike_lane_total.2011ct/100,
            `% Bike lane to streets` = Bike_lane.by.street.2011ct,
            `Green space area (ha)` = area_esp_vert_2011ct,
            `% green space within CT` = pct_esp_vert_2011ct,
            `Tree canopy area (ha)` = area_esp_vert_high_2011ct,
            `% tree canopy within CT` = pct_esp_vert_high_2011ct)

# Extract SES and BEI in 2016
.tab1_2016 <- bei_df_aoi %>%
  transmute(year = "2016",
            quintile = case_when(gentrified_2016_2011 ~ "Gentrified",
                                 TRUE ~ "Not gentrified"),
            zone = zone,
            CT_UID = CT_UID,
            `Bike lane` = Bike_lane_total.2016ct > 0,
            `Bike lane length (hm)` = Bike_lane_total.2016ct/100,
            `% Bike lane to streets` = Bike_lane.by.street.2016ct,
            `Green space area (ha)` = area_esp_vert_2017ct,
            `% green space within CT` = pct_esp_vert_2017ct,
            `Tree canopy area (ha)` = area_esp_vert_high_2017ct,
            `% tree canopy within CT` = pct_esp_vert_high_2017ct)

.tab1_med <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year, quintile) %>%
  summarize(across(where(is.numeric), median, na.rm=TRUE), 
            `% gentrified` = 100 * sum(`gentrified`, na.rm = TRUE) / n(),
            `Bike lane presence` = sum(`Bike lane`))

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

.tab1_p25 <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year, quintile) %>%
  summarize(across(where(is.numeric), quantile, probs=.25, na.rm=TRUE))

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

.tab1_p75 <- bind_rows(.tab1_2011, .tab1_2016) %>%
  group_by(year, quintile) %>%
  summarize(across(where(is.numeric), quantile, probs=.75, na.rm=TRUE))

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

tab1 <- inner_join(.tab1_med, .tab1_p25, by=c("year", "quintile"), suffix = c("", " p25")) %>%
  inner_join(.tab1_p75, by=c("year", "quintile"), suffix = c("", " p75")) %>%
  transmute(year.quintile = paste(year, quintile),
            Population = case_when(year == "2016" ~ "", 
                                   TRUE ~ paste0(round(`Population`,1), " (", round(`Population p25`,1), " - ", round(`Population p75`,1), ")")),
            `Material score` = case_when(year == "2016" ~ "", 
                                         TRUE ~ paste0(round(`Material score`,3), " (", round(`Material score p25`,3), " - ", round(`Material score p75`,3), ")")),
            `% visible minorities` = case_when(is.na(`% visible minorities`) ~ "", 
                                               TRUE ~ paste0(round(`% visible minorities`,1), " (", round(`% visible minorities p25`,1), " - ", round(`% visible minorities p75`,1), ")")),
            `% gentrified CTs` = case_when(year == "2016" ~ "", 
                                           TRUE ~ as.character(round(`% gentrified`, 1))),
            `Bike lane presence (N)` = as.character(`Bike lane presence`),
            `Bike lane length (hm)` = paste0(round(`Bike lane length (hm)`,1), " (", round(`Bike lane length (hm) p25`,1), " - ", round(`Bike lane length (hm) p75`,1), ")"),
            `% Bike lane to streets` = paste0(round(`% Bike lane to streets`,1), " (", round(`% Bike lane to streets p25`,1), " - ", round(`% Bike lane to streets p75`,1), ")"),
            `Green space area (ha)` = paste0(round(`Green space area (ha)`,1), " (", round(`Green space area (ha) p25`,1), " - ", round(`Green space area (ha) p75`,1), ")"),
            `% green space within CT` = paste0(round(`% green space within CT`,1), " (", round(`% green space within CT p25`,1), " - ", round(`% green space within CT p75`,1), ")"),
            `Tree canopy area (ha)` = paste0(round(`Tree canopy area (ha)`,1), " (", round(`Tree canopy area (ha) p25`,1), " - ", round(`Tree canopy area (ha) p75`,1), ")"),
            `% tree canopy within CT` = paste0(round(`% tree canopy within CT`,1), " (", round(`% tree canopy within CT p25`,1), " - ", round(`% tree canopy within CT p75`,1), ")")) %>%
  tibble::column_to_rownames("year.quintile") %>%
  t() %>% 
  as.data.frame()

tab1

9.2 Interactions in UI associations

# extract the legend from one of the plots
legend.row1 <- get_legend(
  # create some space to the left of the legend
  g_mat_bk.hrdl + labs(colour='Material score\n(Quintiles)', fill='Material score\n(Quintiles)') + theme(legend.box.margin = margin(0, 0, 0, 12))
)

prow1.1 <- plot_grid(
  g_mat_bk.hrdl + xlab("% Bike lane to streets 2011") + ylab("% Bike lane to streets 2016") + theme(legend.position = "none"),
  g_mat_gr + xlab("% green space within CT 2011") + ylab("% green space within CT 2016") + theme(legend.position = "none"),
  g_mat_tc + xlab("% tree canopy within CT 2011") + ylab("% tree canopy within CT 2016") + theme(legend.position = "none"),
  align = "vh",
  labels = c("A", "B", "C"),
  hjust = -1,
  nrow = 1
)

prow1 <- plot_grid(prow1.1, legend.row1, rel_widths = c(3, .4))

# extract the legend from one of the plots
legend.row2 <- get_legend(
  # create some space to the left of the legend
  g_vis_bk + labs(colour='Visible minorities\n(Quintiles)', fill='Visible minorities\n(Quintiles)') + theme(legend.box.margin = margin(0, 0, 0, 12))
)

prow2.1 <- plot_grid(
  g_vis_bk.hrdl + xlab("% Bike lane to streets 2011") + ylab("% Bike lane to streets 2016") + theme(legend.position = "none"),
  g_vis_gr + xlab("% green space within CT 2011") + ylab("% green space within CT 2016") + theme(legend.position = "none"),
  g_vis_tc + xlab("% tree canopy within CT 2011") + ylab("% tree canopy within CT 2016") + theme(legend.position = "none"),
  align = "vh",
  labels = c("D", "E", "F"),
  hjust = -1,
  nrow = 1
)

prow2 <- plot_grid(prow2.1, legend.row2, rel_widths = c(3, .4))

# extract the legend from one of the plots
legend.row3 <- get_legend(
  # create some space to the left of the legend
  g_gen_bk + labs(colour='Gentrified', fill='Gentrified') + theme(legend.box.margin = margin(0, 0, 0, 12))
)

prow3.1 <- plot_grid(
  g_gen_bk.hrdl + xlab("% Bike lane to streets 2011") + ylab("% Bike lane to streets 2016") + theme(legend.position = "none"),
  g_gen_gr + xlab("% green space within CT 2011") + ylab("% green space within CT 2016") + theme(legend.position = "none"),
  g_gen_tc + xlab("% tree canopy within CT 2011") + ylab("% tree canopy within CT 2016") + theme(legend.position = "none"),
  align = "vh",
  labels = c("G", "H", "I"),
  hjust = -1,
  nrow = 1
)

prow3 <- plot_grid(prow3.1, legend.row3, rel_widths = c(3, .4))

# output all graphs
plot_grid(prow1, prow2, prow3, ncol=1)

10 Annexes

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] cancensus_0.5.0  DHARMa_0.4.5     stargazer_5.2.2  spaMM_3.10.0    
##  [5] cachem_1.0.6     memoise_2.0.1    spatialreg_1.2-1 spdep_1.2-2     
##  [9] spData_2.0.1     sp_1.4-6         lme4_1.1-28      Matrix_1.4-0    
## [13] cowplot_1.1.1    biscale_0.2.0    openxlsx_4.2.5   RPostgres_1.4.3 
## [17] DBI_1.1.2        ggmap_3.0.0      ggplot2_3.3.5    stars_0.5-5     
## [21] abind_1.4-5      sf_1.0-6         tidyr_1.2.0      dplyr_1.0.8     
## 
## loaded via a namespace (and not attached):
##   [1] minqa_1.2.4           colorspace_2.0-3      rjson_0.2.21         
##   [4] deldir_1.0-6          ellipsis_0.3.2        class_7.3-20         
##   [7] rstudioapi_0.13       proxy_0.4-26          farver_2.1.0         
##  [10] bit64_4.0.5           fansi_1.0.2           lubridate_1.8.0      
##  [13] codetools_0.2-18      splines_4.1.2         doParallel_1.0.17    
##  [16] knitr_1.37            jsonlite_1.8.0        nloptr_2.0.0         
##  [19] png_0.1-7             shiny_1.7.1           compiler_4.1.2       
##  [22] httr_1.4.2            assertthat_0.2.1      fastmap_1.1.0        
##  [25] cli_3.2.0             later_1.3.0           s2_1.0.7             
##  [28] ROI.plugin.glpk_1.0-0 htmltools_0.5.2       tools_4.1.2          
##  [31] coda_0.19-4           gtable_0.3.0          glue_1.6.2           
##  [34] wk_0.6.0              gmodels_2.18.1        Rcpp_1.0.8           
##  [37] slam_0.1-50           jquerylib_0.1.4       raster_3.5-15        
##  [40] vctrs_0.3.8           gdata_2.18.0          nlme_3.1-155         
##  [43] iterators_1.0.14      lwgeom_0.2-8          xfun_0.29            
##  [46] stringr_1.4.0         mime_0.12             lifecycle_1.0.1      
##  [49] gtools_3.9.2          terra_1.5-21          ROI_1.0-0            
##  [52] LearnBayes_2.15.1     MASS_7.3-55           scales_1.1.1         
##  [55] promises_1.2.0.1      hms_1.1.1             parallel_4.1.2       
##  [58] expm_0.999-6          RColorBrewer_1.1-2    yaml_2.3.5           
##  [61] pbapply_1.5-0         sass_0.4.0            stringi_1.7.6        
##  [64] highr_0.9             foreach_1.5.2         gap_1.2.3-1          
##  [67] e1071_1.7-9           boot_1.3-28           zip_2.2.0            
##  [70] qgam_1.3.4            RgoogleMaps_1.4.5.3   rlang_1.0.1          
##  [73] pkgconfig_2.0.3       bitops_1.0-7          evaluate_0.15        
##  [76] lattice_0.20-45       purrr_0.3.4           labeling_0.4.2       
##  [79] bit_4.0.4             tidyselect_1.1.2      plyr_1.8.6           
##  [82] magrittr_2.0.2        geojsonsf_2.0.1       R6_2.5.1             
##  [85] generics_0.1.2        mgcv_1.8-39           pillar_1.7.0         
##  [88] withr_2.4.3           units_0.8-0           tibble_3.1.6         
##  [91] crayon_1.5.0          Rglpk_0.6-4           KernSmooth_2.23-20   
##  [94] utf8_1.2.2            rmarkdown_2.11        jpeg_0.1-9           
##  [97] grid_4.1.2            blob_1.2.2            digest_0.6.29        
## [100] classInt_0.4-3        xtable_1.8-4          httpuv_1.6.5         
## [103] numDeriv_2016.8-1.1   munsell_0.5.0         registry_0.5-1       
## [106] bslib_0.3.1